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

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 = SimplifyPWToUnitStep@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]] &;

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!

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]] &;

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.

• Code you posted has several typos, not only with q2. Please, clean up all of them from your new post. Dec 26, 2022 at 5:40
• Hello! Thanks for your help and I've looked through the code again and edited all the typos I am able to identify... But honestly I've been quite bad at this, if you seen anything else please feel free to notify me as well!
– FLP
Dec 26, 2022 at 11:06
• Also the definition of pdetoode is not included directly in the code, so if you are running it you would have to include it separately - not sure if this would help as well
– FLP
Dec 26, 2022 at 11:07
• Actually numerical model with rb = 100; points = 200; difforder = 2; is very rough. I can recommend to map solution on the unit interval with lb=10^-3, $lb\le r\le 1$ and normalize h1, h2 as h1/100, h2/100. Then take points=30` in the first run. Dec 27, 2022 at 3:49
• It is not clear why in your plot $r<0$ while it should be $10^{-3}\le r \le 1$? Dec 28, 2022 at 6:00