I'm trying to numerically solve a system of PDEs, under some analytically derived initial conditions and matching boundary conditions - i.e. just forcing the BCs to be the ICs evaluated at the endpoints (the code is longish, so I'll include it at the bottom).
However, the IC for one of the variables t
, is quite explosive on its one end point:
which, I believe is the reason why I've been getting the inconsistent ICs/BCs warning (NDSolve::ibcinc) - since they just can't be inconsistent, and I've been able to get rid of the warning by adjusting the spatial discretization (as I saw recommended on this site).
Assuming that I'm right about that, I don't know how to go about knowing how best to tinker with the discretization. So far, I've been able to get rid of that inconsistency warning by taking the number of points used by "TensorProductGrid"
to ~10^5 - which seems absurd, and also means that I'm not able to temporally evolve the solution for longer than ~10^-14 without it crashing either my kernel or my whole laptop. I've also been able to get rid of the warning by taking the "TensorProductGrid"
number of points to ~833, but with the "DifferenceOrder"
set to "Pseudospectral"
(which I've seen advised, but don't understand since it seems like this should work only with periodic BCs). But, I'm still not able to fully evolve the solution temporally, and it starts introducing stiffness warnings (NDSolve::ndsz) that I haven't been able to deal with.
Has anyone got any idea as to what I should be doing? I haven't been able to find any fitting advice on here yet.
Here's the code:
T = 1/(π uh);
L = 1;
rh = 1; rc = 10^3;
rs = rh + 10^-3;
uc = 1/rc;
uh = 1/rh;
l0 = rc - rs;
x0 = rs;
a = 2.6;
b = 2;
ainc = 0.2;
τmin = 0;
τmax = 1;
f[u_, uhpar_] := 1 - (u/uhpar)^2
Σ[u_, ucpar_, uhpar_, ap_, bp_] :=
((1 - uhpar/u)/(1 - uhpar/ucpar))^ap (u/ucpar)^bp;
pde[ucpar_, uhpar_, ap_, bp_] :=
{D[(1/Σ[u[τ, σ], ucpar, uhpar, ap, bp])((f[u[τ, σ], uhpar]*L^2)/u[τ, σ]^2)D[t[τ, σ], τ], τ] - D[Σ[u[τ, σ], ucpar, uhpar, ap, bp]((f[u[τ, σ], uhpar] L^2)/u[τ, σ]^2)D[t[τ, σ], σ], σ] == 0,
D[(1/Σ[u[τ, σ], ucpar, uhpar, ap, bp]) (L^2/u[τ, σ]^2)D[x[τ, σ], τ], τ] - D[Σ[u[τ, σ], ucpar, uhpar, ap, bp] (L^2/u[τ, σ]^2)D[x[τ, σ], σ], σ] == 0,
D[(1/Σ[u[τ, σ], ucpar, uhpar, ap, bp]) (L^2/u[τ, σ]^2)(1/f[u[τ, σ], uhpar])D[u[τ, σ], τ], τ] - D[Σ[u[τ, σ], ucpar, uhpar, ap, bp](L^2/u[τ, σ]^2)(1/f[u[τ, σ], uhpar])D[u[τ, σ], σ], σ] ==
(-(1/2))*(1/Σ[u[τ, σ], ucpar, uhpar, ap, bp])((D[t[τ, σ], τ]^2 - Σ[u[τ, σ], ucpar, uhpar, ap, bp]^2 D[t[τ, σ], σ]^2)D[(f[u[τ, σ], uhpar]L^2)/u[τ, σ]^2, u[τ, σ]] - (D[x[τ, σ], τ]^2 - Σ[u[τ, σ], ucpar, uhpar, ap, bp]^2 D[x[τ, σ], σ]^2) D[L^2/u[τ, σ]^2, u[τ, σ]] - (D[u[τ, σ], τ]^2 - Σ[u[τ, σ], ucpar, uhpar, ap, bp]^2*D[u[τ, σ], σ]^2) D[(L^2/u[τ, σ]^2)(1/f[u[τ, σ], uhpar]), u[τ, σ]])}
rsStar = L^2/rh ArcCoth[-(rs/rh)];
σf = L^2/rh ArcCoth[-((rs + l0)/rh)] - rsStar;
rAnalytical[σ_] := -rh Coth[(rh(rsStar + σ))/L^2];
tInit[σ_] := σ - Coth[σ - ArcCoth[1001/1000]];
dtInit[σ_] := Coth[σ - ArcCoth[1001/1000]]^2;
newBegBC =
{t[τ, σ] == tInit[0], D[x[τ, σ], σ] == 0, u[τ, σ] == rs^(-1)} /. σ -> 0;
newEndBC =
{D[t[τ, σ], σ] == dtInit[σf], D[x[τ, σ], σ] == 0,
D[u[τ, σ], σ] == (D[rAnalytical[tσ]^(-1), tσ] /. tσ -> σf)
newIC =
{t[tp, σ] == tInit[σ],
x[tp, σ] == 0,
u[tp, σ] == rAnalytical[σ]^(-1),
(D[t[τ, σ], τ] /. τ -> tp) == 0,
(D[x[τ, σ], τ] /. τ -> tp) == 0,
(D[u[τ, σ], τ] /. τ -> tp) == 0} /. tp -> τmin;
(* attempt 1 *)
ngrd = 9 10^5;
NDSolve[
Flatten[{pde[uc, uh, a, b], newBegBC, newEndBC, newIC}],
{t, x, u}, {τ, 0, 10^-14}, {σ, 0, σf},
PrecisionGoal -> 3, AccuracyGoal -> 4,
Method ->
{"MethodOfLines",
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1},
"SpatialDiscretization" ->
{"TensorProductGrid",
"MaxPoints" -> ngrd,
"MinPoints" -> ngrd,
"DifferenceOrder" -> 4}}]
(* attempt 2 *)
ngrd = 833;
NDSolveValue[
Flatten[{pde[uc, uh, a, b], newBegBC, newEndBC, newIC}],
{t, x, u}, {τ, min, 10^-12}, {σ, 0, σf},
PrecisionGoal -> 2, AccuracyGoal -> 2,
Method ->
{"MethodOfLines",
"SpatialDiscretization" ->
{"TensorProductGrid",
"MaxPoints" -> ngrd,
"MinPoints" -> ngrd,
"DifferenceOrder" -> "Pseudospectral"}}]