# 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. – Feyre Jun 14 '16 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? – dearN Jun 14 '16 at 15:26
• How about just increasing precision? PrecisionGoal -> 10 gives an error of -0.470581. – Feyre Jun 14 '16 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? – dearN Jun 14 '16 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. – Feyre Jun 14 '16 at 16:32

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.