Speeding up NDSolve to reasonable speeds to solve a coupled PDE system

Question

Problem Statement

I am planning to solve a PDE system which consists of a fluid droplet spreading on a non-Newtonian substrate. The system consists of the following equations:

$$\frac{\partial p_1}{\partial r}=\kappa (\frac{\partial u_1}{\partial z})^n\rightarrow u_1=\frac{(\frac{z}{\kappa}\frac{\partial p_1}{\partial r}+c_1)^{1+\frac{1}{n}}}{\frac{1}{\kappa}\frac{\partial p_1}{\partial r}(1+\frac{1}{n})}+c_2$$ $$\frac{\partial p_2}{\partial r}=\mu \frac{\partial^2 u_2}{\partial z^2} \rightarrow u_2=\frac{z^2}{2\mu}\frac{\partial p_2}{\partial r}+a_1z+a_2$$

$$\frac{\partial h_1}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\int_{0}^{h_1} r u_1 \,dz$$ $$\frac{\partial h_2}{\partial t}=\frac{1}{r}\frac{\partial}{\partial r}\int_{h_1}^{h_2} r u_2 \,dz$$

The constants $c_1,c_2,a_1,a_2$ can be solved by applying the following boundary conditions: $$p_1=-\gamma_1 (\frac{\partial^2 h_1}{\partial r^2}+\frac{1}{r}\frac{\partial h_1}{\partial r})-\gamma_2 (\frac{\partial^2 h_2}{\partial r^2}+\frac{1}{r}\frac{\partial h_2}{\partial r})$$ $$p_2=-\gamma_2 (\frac{\partial^2 h_2}{\partial r^2}+\frac{1}{r}\frac{\partial h_2}{\partial r})-A(\frac{h_f}{h})^3$$

At $z=h_1,\ \kappa (\frac{\partial u_1}{\partial z})^n=\mu \frac{\partial^2 u_2}{\partial z^2}$

At $z=h_2,\ \mu \frac{\partial^2 u_2}{\partial z^2}=0$

Ive managed to construct a code to solve for the combined system of equations using pdetoode. It solves at a fine steady rate, but as of now its solving speed is really really low at around $10^{-4}$ simulated time per second.

As I'll need to solve the equations to a time of $10 ^3$, at the current rate I'll have to wait for 115 days, which is nearly half a year D: Clearly that is not sustainable at all.

I'm wondering if anyone here knows alternative methods that can speed up the code to a manageable rate. Thanks everyone to all the help offered!

Code

(*Defining constants*)
(It is the current limitation of the code that n can only be either 

odd or 1/odd)
n = 1/3(shear thinning power);
[Kappa] = 1(stress constant);
[Mu] = 1(viscosity);
G1 = 0.01(surface tension of bottom fluid);
G2 = 1(surface tension of top fluid);
a = 0.01(disjoined pressure constant);
hf = 0.01(height of disjoined pressure);
[Epsilon] = 1*10^-9;
[CapitalPi][h_] := a (hf/h)^3(*(1-(hf/h)^3));
(Inputting Equations)
unitStepExpand = Simplify`PWToUnitStep@PiecewiseExpand@# &;
With[{p1 = p1[r, t], p2 = p2[r, t], h1 = h1[r, t], h2 = h2[r, t], 
  q1 = q1[r, t], q2 = q2[r, t], c1 = c1[r, t], c2 = c2[r, t], 
  a1 = a1[r, t], a2 = a2[r, t]},
(Defining pressure)
 Eqnp1 = p1 == -G2 (D[h2, r, r] + 1/r D[h2, r]) - 
    G1 (D[h1, r, r] + 1/r D[h1, r]);
 Eqnp2 = p2 == -G2 (D[h2, r, r] + 1/r D[h2, r]) - [CapitalPi][h2];
(Defining the velocity constants)
 c1 = 1/[Kappa] ((h1 - h2) D[p2, r] - h1 D[p1, r]);
 c2 = -(c1^(1 + 1/n)/(1/[Kappa] D[p1, r] (1 + 1/n)));
 a1 = -h2/[Mu] D[p2, r];
 a2 = (h1/[Kappa] D[p1, r] + c1)^(1 + 1/n)/(
   1/[Kappa] D[p1, r] (1 + 1/n)) + c2 - h1^2/(2 [Mu]) D[p2, r] - 
   a1 h1;
(Defining velocity and flux)
 u1 = (z/[Kappa] D[p1, r] + c1)^(1 + 1/n)/(
   1/[Kappa] D[p1, r] (1 + 1/n)) + c2;
 u2 = z^2/(2 [Mu]) D[p2, r] + a1 z + a2;
Eqnq1 = q1 == 
   Simplify[
    r Integrate[u1, {z, 0, h1}] /. 
     D[p1, r] -> 
      unitStepExpand@
       If[Sqrt[D[p1, r]^2] < [Epsilon], [Epsilon], D[p1, r]]];
 Eqnq2 = q2 == 
   Simplify[
    r Integrate[u2, {z, h1, h2}] /. 
     D[p1, r] -> 
      unitStepExpand@
       If[Sqrt[D[p1, r]^2] < [Epsilon], [Epsilon], D[p1, r]]];
(The unitstepexpand functions are needed to prevent q1 and q2 from being evaluated as complexinfinity at the beginning of the solver, since the pressure and their respective derivatives evaluates to 0)
(Defining final equations)
Eqnh1 = D[h1, t] + D[q1, r] == 0;
 Eqnh2 = D[h2, t] + D[q1, r] + D[q2, r] == 0;
]
(Defining Initial and Boundary Conditions)
hbath = 100;
IC = {h1[r, 0] == hbath - 10 E^-(r/10)^2, h2[r, 0] == hf + hbath}
BC = {{D[h1[r, t], r] == 0, D[h1[r, t], r, r, r] == 0, 
    D[h1[x, t], x] == 0, D[h1[x, t], x, x, x] == 0, 
    D[h2[r, t], r] == 0, D[h2[r, t], r, r, r] == 0, 
    D[h2[x, t], x] == 0, D[h2[x, t], x, x, x] == 
          0} /. {r -> lb, x -> rb}}
(Generation of grid)
lb = 1/1000000; rb = 100; points = 200; difforder = 2;
grid = Array[# &, points, {lb, rb}];
(Discretisation)
ptoofunc = 
  pdetoode[{p1, p2, h1, h2, q1, q2}[r, t], t, grid, difforder];
removeredundant = #[[3 ;; -3]] &;
odeadd = Map[
   ptoofunc, {Eqnp1, Eqnp2, Eqnq1, 
    Eqnq2}, {2}];
ode1 = Block[{p1, p2, h1, h2, q1, q2}, Set @@@ odeadd; 
    ptoofunc@Eqnh1] // removeredundant;
ode2 = Block[{p1, p2, h1, h2, q1, q2}, Set @@@ odeadd; 
    ptoofunc@Eqnh2] // removeredundant;
ode = {ode1, ode2};
odIC = ptoofunc@IC;
With[{sf = 1}, odBC = diffbc[t, sf]@BC // ptoofunc];
(Solving)
var = Outer[#[#2] &, {h1, h2}, grid];
Monitor[sollst = 
   NDSolveValue[{ode, odIC, odBC}, var, {t, 0, 1000}, 
    Method -> {"EquationSimplification" -> "MassMatrix"}, 
    EvaluationMonitor :> (time = t)], time];

Edits

I've tried all optimisation methods I knew, such as adjusting AccuracyGoal, PrecisionGoal to 1, and even reduced the number of points solved to 20, and yet the solution does not solve at any reasonable rate at all, its still stuck between $10^{-4}$ to $10^{-6}$ simulated time per second... Ive never seen a code so stubbornly slow...

Edits 2

Edited all the typos I was able to find, if there is anything else please notify me!

Answer 1

Apologies but apparently the error lies in a typo in the equations:

 (*Defining final equations*)
Eqnh1 = D[h1, t] + D[q1, r] == 0;
 Eqnh2 = D[h2 - h1, t] + D[q1, r] == 0;

Comparing the 2nd final equation with the supposed one, we notice that Eqnh2 is of the wrong form and should be

 (*Defining final equations*)
Eqnh1 = D[h1, t] + D[q1, r] == 0;
 Eqnh2 = D[h2 - h1, t] + D[q2, r] == 0;

So sorry for everyone who have tried and help out with the problem... Personally I didn't expect this to be the cause of the error since the code manages to solve fine, and I was never able to observe any numerical results due to the short time frame of the code.

The code of course still faces other issues.

Plotting the pressure of the non-Newtonian layer, we notice that the pressure profile has a sudden spike as $r\rightarrow 0 $.

At later times, the error further propagate down the solution.

Since the solution is supposed to be axisymmetric, there is no physical reason why the pressure would spike at $r\rightarrow 0 $. I think it is likely cause by my choice of boundary conditions

BC = {{D[h1[r, t], r] == 0, D[h1[r, t], r, r, r] == 0, 
    D[h1[x, t], x] == 0, D[h1[x, t], x, x, x] == 0, 
    D[h2[r, t], r] == 0, D[h2[r, t], r, r, r] == 0, 
    D[h2[x, t], x] == 0, D[h2[x, t], x, x, x] == 
          0} /. {r -> lb, x -> rb}}

Since the left boundary is not exactly 0, enforcing that the first and third derivatives = 0 may have caused the unphysical spike in pressure.

Ive therefore changed the boundary condition of the code to Periodic, enforced by the True option in pdetoode, and solved the problem in full space instead of half space:

(*Generation of grid*)lb =-100; rb = 100; points = 200; difforder = 2;
grid = Array[# &, points, {lb, rb}];
(Discretisation)
ptoofunc = 
  pdetoode[{p1, p2, h1, h2, q1, q2}[r, t], t, grid, difforder, True];
removeredundant = #[[3 ;; -3]] &;
odeadd = Map[
   ptoofunc, {Eqnp1, Eqnp2, Eqnq1, 
    Eqnq2}, {2}];
ode1 = Block[{p1, p2, h1, h2, q1, q2}, Set @@@ odeadd; 
    ptoofunc@Eqnh1] (// removeredundant);
ode2 = Block[{p1, p2, h1, h2, q1, q2}, Set @@@ odeadd; 
    ptoofunc@Eqnh2] (// removeredundant);
ode = {ode1, ode2};
odIC = ptoofunc@IC;
(With[{sf = 1}, odBC = diffbc[t, sf]@BC // ptoofunc];)
(Solving)
var = Outer[#[#2] &, {h1, h2}, grid];
Monitor[sollst = 
   NDSolveValue[{ode, odIC (odBC)}, var, {t, 0, 1000}, 
    Method -> {"EquationSimplification" -> "MassMatrix"}, 
    EvaluationMonitor :> (time = t)], time];

While it solves fine it stiffens rapidly at time scale of $10^{-7}$, however since it is a separate numerical issue I would leave it for another post.

Edit

Since theres interest in the problem I have specified the exact changes I have made that have led to the current numerical instability, hope this would help!

Edit 2

I have normalised the values of h1, h2 and ran the solution again at the advice of Alex Trounev. Here are the results:

Again, the pressure singularity at the centre seems to be the culprit:

However, even after editing out the pressure to remove the pressure singularity:

(*Defining pressure*)
Eqnp1 = p1 == -G2 (D[h2, r, r]) - 
    G1 (D[h1, r, r](*+(1/r)D[h1,r]*));
Eqnp2 = p2 == -G2 (D[h2, r, r]) - \[CapitalPi][h2];

The pressure field continue to explode at the left hand side - it seems that for some reason the left hand side of the plot is surprisingly unstable.