# NDSolve struggling with tricky boundary conditions

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[τ, σ]])}

σ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, 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"}}]


Reasonably accurate discretization of tInit[σ], the initial condition for t, requires that it not change much from one grid point to the next. The rate of change of tInit[σ] is its derivative,

slope = D[tInit[σ], σ]
(* 1 + Csch[σ - ArcCoth[1001/1000]]^2 *)


which varies by six orders of magnitude across the domain of integration. Hence, limiting change per grid point to, say, 10%, would require 10^7 grid points on a uniform mesh. It is, however, not necessary to use a uniform mesh. Instead, use

grid = Join[{0}, Table[σ /. FindRoot[tInit[σ] == i 10^-4 (tInit[σf] - tInit) + tInit,
{σ, 0, σf}], {i, 9999}], {σf // N}];


which limits the variation of tInit[σ] from grid point to grid point,

MinMax@Differences[Table[tInit[σ], {σ, grid}]]
{0.10028, 0.10028}


as desired. Directing NDSolve to use this grid is accomplished by

s = Flatten@NDSolve[{pde[uc, uh, a, b], newBegBC, newEndBC, newIC}, {t, x, u},
{τ, 0, 10^-9}, {σ, 0, σf}, Method -> {"MethodOfLines", "SpatialDiscretization" ->
{"TensorProductGrid", "Coordinates" -> {grid}}, Method -> "StiffnessSwitching"}]


which permits the computation to proceed five orders of magnitude further than in the question. (The error message described in the question also goes away.) Proceeding further is limited by memory usage or stiffness, depending on other options used by NDSolve. Of course, this is not nearly far enough, as can be seen by comparing solutions at τ == 0 and τ == 10^-9.

LogPlot[Evaluate[{t[0, σ], t[10^-9, σ]} /. s], {σ, 0, σf}, PlotRange -> All,
ImageSize -> Large, AxesLabel -> {σ, t}, LabelStyle -> Directive[12, Black, Bold]] The two curves are indistinguishable. Perhaps, stiffness can be circumvented by some other approach. It seems unlikely, however, that spatial discretization can be improved further.