# Speeding up WhenEvent + NDSolve (Performance tuning and handling instabilities)

I am solving a non-linear partial differential equation using NDSolve. The equation models a liquid film on an isothermal solid substrate. As a result of various fluid physics mechanisms such as gravity, surface tension, thermocapillarity, this liquid film contorts, wrinkles and eventually "ruptures".

I am trying to use WhenEvent to StopIntegration in NDSolve when the film thickness, u is below a certain value (0.0001, in this case). This is meant to be accomplished by the following WhenEvent statement:

 WhenEvent[Min[ParallelTable[u[t, x], {x, xMin, xMax}]] < 0.0001,
end = t; "StopIntegration"]


The Mathematica code for this is:

Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
j = 2;
k = 0.0677;
{xMin, xMax} = {0, j*(Pi/k)};
TMax = 200000;
Bi = 1.; M = 35.1; Pr = 7.02; Ga = 0.; mult = 1.; Pm = 0.025; S = 100.;
AbsoluteTiming[
Quiet[hSol =
u /. NDSolve[{D[u[t, x],
t] == (-S)*D[u[t, x]^3*D[u[t, x], x, x, x], x] + (Ga/3)*
D[u[t, x]^3*D[u[t, x], x], x] -
(Bi*(M/Pr))*
D[(u[t, x]/(1 + Bi*u[t, x]))^2*D[u[t, x], x], x] +
mult*Pm*(D[u[t, x], x, x] - Ga*u[t, x]),
u[0, x] == 1 + 0.1*Cos[k*x], u[t, xMin] == u[t, xMax],
Derivative[0, 1]*u[t, xMin] == Derivative[0, 1]*u[t, xMax],
Derivative[0, 2]*u[t, xMin] == Derivative[0, 2]*u[t, xMax],
Derivative[0, 3]*u[t, xMin] == Derivative[0, 3]*u[t, xMax],
WhenEvent[
Min[ParallelTable[u[t, x], {x, xMin, xMax}]] < 0.0001,
end = t; "StopIntegration"]},
u, {t, 0, TMax}, {x, xMin, xMax}, MaxSteps -> 100000,
PrecisionGoal -> Automatic,
Method -> {"MethodOfLines", "Method" -> "LSODA",
"TemporalVariable" -> t, "SpatialDiscretization" ->
{"TensorProductGrid", "MinPoints" -> 800,
"MaxPoints" -> 1200, "DifferenceOrder" -> 5}}][[1]]]]


NDSolve takes ~1 minute on my laptop which is Intel Core i7-2670QM, 2.20 GHz x 8 with 8 GB memory.

However, I notice that when WhenEvent ends the simulation, I have a huge error in mass conservation:

 hGrid = InterpolatingFunctionGrid[hSol];
tRup = InterpolatingFunctionDomain[hSol][[1]][[2]];
mass = NIntegrate[1 + 0.1*Cos[k*x], {x, xMin, xMax},
Method -> {Automatic, "SymbolicProcessing" -> 0}]

errorInMassCon=((NIntegrate[hSol[tRup, x], {x, xMin, xMax},
Method -> {Automatic, "SymbolicProcessing" -> 0}]) - mass)*100/
mass (*14%*)


In the above simulation, when I set the multiplier, mult=0. instead of mult=1., it makes a big difference to the solution time (complete in ~12 seconds) and proper mass conservation. The reduction in solution time is expected as setting mult=0. drops out the following terms: mult*Pm*(D[u[t, x], x, x] - Ga*u[t, x]).

mass = NIntegrate[1 + 0.1*Cos[k*x], {x, xMin, xMax},
Method -> {Automatic, "SymbolicProcessing" -> 0}];

errorInMassCon=((NIntegrate[hSol[tRup, x], {x, xMin, xMax},
Method -> {Automatic, "SymbolicProcessing" -> 0},
MaxRecursion -> 15]) - mass)*100/mass (*~1%*)


Is there a way to tune the performance so that I can speed up the simulation and avoid this instability (mass conservation failure), if this is a numerical artifact?

Here is what the wrinkled surface of the film looks like as a result of an interplay between various mechanisms (to those interested :)):

A schematic of various mechanisms that affect liquid film dynamics:

• Hmm, in both cases I get 69.1977 for errorInMassCon. It takes my rig approx 17 secs to run the original one with 16g ram, huge performance boost which isn't expensive, I've found a fast usb stick with readyboost helps with my previous 6g ram laptop. Jun 14, 2016 at 14:51
• @Feyre 1) I tried this on a 24gig machine, the issue with mass conservation error still exits. There must be some sort of "numerical viscosity" that can be included as a tuning parameter for a stiff problem of this fashion? Jun 14, 2016 at 15:26
• How about just increasing precision? PrecisionGoal -> 10 gives an error of -0.470581. Jun 14, 2016 at 15:56
• @Feyre That (PrecisionGoal->10) seems to have helped with mult=1.. However, whats wrong with PrecisionGoal->Automatic? How does one tweak these parameters? Trial-and-error? Jun 14, 2016 at 16:07
• Automatic Gives the precision goal of half of MachinePrecision, that doesn't mean it's optimized for your calculation. I just increased it incrementally. When you have high error, precision is the first you look at right? I honestly don't know why that term is specifically giving such problems though. It's possible that adding 2nd derivative limits the amount of solution functions, change the order of differentiation there and the required time plummets. Jun 14, 2016 at 16:32

## 1 Answer

PrecisionGoal -> 10 gives an error of -0.470581 in the weight.

Automatic Gives the precision goal of half of MachinePrecision, that doesn't mean it's optimized for your calculation. When you have high error, precision is the first you look at, so incrementally increasing precision can help.

It's possible that adding 2nd derivative limits the amount of solution functions, change the order of differentiation there and the required time plummets. This means that it's the complexity of the ODE that makes the calculations take so long., and the NDSolve itself probably can't be further optimized.