Difficulty in solving a differential equation with custom solver for a different initial condition

Question

This is a follow-up question to this previous post by @FLP, in which an interesting system of equations was solved with the useful pdetoode developed by @xzczd. I have tried to solve this problem with a different initial condition (IC). And I appreciated the help and effort from @xzczd and @FLP. The following code is a combination of the original post and the answer. Due to the symmetry, the differential equation for h[r,t] is to be solved from $r=0$ to a sufficiently large value $r>>1$. Here, the left bound of $r$ is taken as a small non-zero value due to some numerical reason, and the right bound is taken as $50$.

My question:

In the original post, a system about $h[r,t]$ was solved with an IC of the exponential function IC = {h[r, t] == 20 E^-(r/20)^2 + hf} /. t -> 0;. For some weird reason if you change the amplitude of the IC, say, IC = {h[r, t] == 10 E^-(r/10)^2 + hf} /. t -> 0;, it would stuck at a small time (relative to a given Tend) and return an error: step size is effectively zero; singularity or stiff system suspected.

I also tried a parabolic IC defined by If:

IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0

but the error of a singularity remained. Noted that in this case PiecewiseExpand transforms If[...] to an equivalent Piecewise[...]; then SimplifyPWToUnitStep transforms Piecewise[...] to an equivalent combination of UnitStep[...]. This is a necessary transformation since pdetoode assumes the functions in the equations, ICs, and BCs are all Listable.

For the piecewise IC, I guess the issue may be caused by the discontinuous derivative at $r=1$. But I cannot understand why the code fails with the IC only changed in its amplitude: h[r, t] == 10 E^-(r/10)^2 + hf.

An even more strange issue is that if I change the side of the domain with, say rb=100, a singularity occurs again.

The code:

The code for pdetoode neglected here, please find in the abovementioned link. Here, Int[r,t] is an integral function, defined as $$Int(r,t)=\int_0^\infty K(r, r^\prime)\frac{\partial}{\partial r^\prime }\left(\frac{h_f}{h} \right)^3 \mathrm{d}r^\prime, $$ where $K(r,r^\prime)=\frac{2}{\pi}r[\mathbf{K}(r^\prime/r)-\mathbf{E}(r^\prime/r)]$ for $r^\prime<r$

or $K(r,r^\prime)=\frac{2}{\pi}r^\prime[\mathbf{K}(r/r^\prime)-\mathbf{E}(r/r^\prime)]$ for $r^\prime>r$, in which $\mathbf{K}$ and $\mathbf{E}$ are the complete elliptic integrals.

In order to avoid the singularity of $K(r,r^\prime)$ at $r=r^\prime$, $r$ takes the values at the grid point lb+i*step and $r^\prime$ takes the values at the staggered grid lb+(i+0.5)*step.

(*constants*)
a = 0.001; hf = 0.01; \[Mu] = 1; \[Beta] = 1;
(*Equations*)
With[{h = h[r, t], P = P[r, t], Q = Q[r, t], J = J[r, t]},
EqnP = P == D[h, r, r] + D[h, r]/r + a^2/h^3;
EqnQ = Q == 1/3*h^3*r*D[P, r];
EqnC = D[h, t] + (\[Mu] D[Q, r])/r == -J;
EqnJ = J == (\[Beta] D[Int[r, t], r])/r;]
(Computation Domain)
lb = 1/100; rb = 50;
(unitStepExpand is a function that transforms If[[Ellipsis]] to mathematically equivalent combination of UnitStep[[Ellipsis]], useful in the computation of the following integral)
unitStepExpand = Simplify`PWToUnitStep[PiecewiseExpand[#]] &;
(initial and boundary conditions)
IC = {h[r, t] == 20 E^-(r/20)^2 + hf} /. t -> 0;
(IC = {h[r, t] == 10 E^-(r/10)^2 + hf} /. t -> 0;)
(IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5(1 - r^2) + hf, hf]} /. t -> 0;*)
BC = {{D[h[r, t], r] == 0, D[h[r, t], r, r, r] == 0, D[h[x, t], x] == 0, D[h[x, t], x, x, x] == 0} /. {r -> lb, x -> rb}};
(Discretization setting)
difforder = 2;
points = 200;
grid = Array[# &, points, {lb, rb}];
hgrid = Map[h[#][t] &, grid];
hdgrid = fdd[Derivative[1], grid, hgrid];
(the Integral Function)
Kernel[r_, i_] := (2r)/[Pi]unitStepExpand@If[r > i, r (EllipticK[i/r] - EllipticE[i/r]), i (EllipticK[r/i] - EllipticE[r/i])];
Integrand = D[(hf/h[r, t])^3, r] /.Thread[{D[h[r, t], r], h[r, t]} -> {hdgrid, hgrid}];
step = (rb - lb)/points;
Int[r_][t_] := Total@Table[stepKernel[r, lb + (i + 0.5)step]*0.5 (Integrand[[i]] + Integrand[[i + 1]]), {i, 1, points - 1}];
(Discretising the rest of the system)
ptoofunc = pdetoode[{h[r, t], P[r, t], Q[r, t], J[r, t], Int[r, t]}, t, grid, difforder];
removeredundant = #[[3 ;; -3]] &;
odeadd = Map[ptoofunc, {EqnP, EqnQ, EqnJ}, {2}];
ode = Block[{P, Q, J}, Set @@@ odeadd; ptoofunc@EqnC] // removeredundant;
odIC = ptoofunc@IC;
With[{sf = 1}, odBC = diffbc[t, sf]@BC // ptoofunc];
(Solving the System)
var = Map[h, grid];
showStatus[status_] := LinkWrite[$ParentLink,SetNotebookStatusLine[FrontEnd`EvaluationNotebook[],ToString[status]]];
clearStatus[] := showStatus[""];
clearStatus[]
monitor[t_] := EvaluationMonitor :> showStatus["t = " <> ToString[CForm[t]]]
Tend = 200;
sollst = NDSolveValue[{ode, odIC, odBC}, var, {t, 0, Tend}, monitor[t], Method -> {EquationSimplification -> MassMatrix}]; // AbsoluteTiming
(postprocess)
solfunc = rebuild[sollst, grid];
Manipulate[Plot[solfunc[t, r], {r, lb, rb}, PlotRange ->All], {t, 0, Tend}]

Some notes:

1. This code works well with the original IC, but when you change the amplitude of the IC, say, from $20$ to $10$, a singularity will be encountered. Increasing the point number and order of FD has little improvement.

2. I am wondering if a nonuniform mesh (or adaptive mesh) could be implemented in the framework of pdetoode such that a refined resolution is bunched towards the region where the change of slope of h[r,t] is large based on physical reason?

3. As mentioned by @FLP, the original code has some additional problems as well. I'm not sure if the treatment of IC will be one of the potential issues. So I post this question here to draw attention.

Please give me some suggestions. Thank you in advance!

Update:

1. With the suggestion of xzczd, the code can run with IC = {h[r, t] == 10 E^-(r/20)^2 + hf} /. t -> 0;, which only changes in its amplitude, up to T=502 approximately, and then singularity occurs at the radius position of the blob, where the slope changes most rapidly.

2. By changing the width of the initial blob, I found smaller the initial width earlier the singularity occurs. For example, using IC = {h[r, t] == unitStepExpand@If[r >= lb && r < 1, 3/5*(1 - r^2) + hf, hf]} /. t -> 0, the singularity occurs at t=0.000024 approximately. Increasing points improves things a little.

Answer 1

If "the IC only changed in its amplitude", changing sf = 1 to sf = 10 resolves the problem.