# Improving NDSolve speed for heavily stiff problems

Having looked around the intergoogles and Mathematica.SE, I thought I'd pose a question with a minimum working example.

Here is the situation I am trying to improve:

1. I am solving a 4th order non linear PDE with NDSolve.
2. It is stiff and I use a stiff solver such as BDF or LSODA.
3. On occassion, I have no choice but to increase the MaxStepFraction to uncomfortable levels.
4. As a result, the code runs longer than usual (made worse by the fact that it is a stiff equation to begin with)

Is there any way I could improve NDSolve performance/speed?

Here is my minimum example:

\$HistoryLength = 0;
Needs["VectorAnalysis"]
Needs["DifferentialEquationsInterpolatingFunctionAnatomy"];
Clear[Eq0, EvapThickFilm, h, Bo, \[Epsilon], K1, \[Delta], Bi, m, r]
Eq0[h_, {Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_}] := \!$$\*SubscriptBox[\(\[PartialD]$$, $$t$$]h\) +
m (h/(K1 + Bi h))^2 Grad[h]] + \[Epsilon]/(
Bi h + K1) + (r) D[D[(h^2/(K1 + Bi h)), x] h^3, x] == 0;
SetCoordinates[Cartesian[x, y, z]];
EvapThickFilm[Bo_, \[Epsilon]_, K1_, \[Delta]_, Bi_, m_, r_] :=
Eq0[h[x, y, t], {Bo, \[Epsilon], K1, \[Delta], Bi, m, r}];
TraditionalForm[EvapThickFilm[Bo, \[Epsilon], K1, \[Delta], Bi, m, r]];
L = 2*92.389;

TMax = 3100*100;
Off[NDSolve::mxsst];
Clear[Kvar];
Kvar[t_] :=  Piecewise[{{1, t <= 1}, {2, t > 1}}]
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol = h /. NDSolve[{
(*Bo,\[Epsilon],K1,\[Delta],Bi,m,r*)

EvapThickFilm[0.003, 0, 1, 0, 1, 0.025, 0],
h[0, y, t] == h[L, y, t],
h[x, 0, t] == h[x, L, t],
(*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)

h[x, y, 0] ==
1 + (-0.25 Cos[2 \[Pi] x/L] - 0.25 Sin[2 \[Pi] x/L]) Cos[
2 \[Pi] y/L]
},
h,
{x, 0, L},
{y, 0, L},
{t, 0, TMax},
Method -> {"BDF", "MaxDifferenceOrder" -> 1},
MaxStepFraction -> 1/50
][] // AbsoluteTiming


A BDF limited to Order 1 needs about 41 seconds to solve the equation until ****failure**** while the LSODA allowed up to order 12 does a fantastic job of cutting it down to 18 seconds.

Now when I increase the MaxStepFraction, I obviously increase the grid density. I am currently running a case that has several thousand grid points and has taken 24+ HOURS, yes hours and hasn't given me a solution yet. This was expected as I have run cases that took about 3-4 hours before with a coarser grid and do hog the ram (they take up about ~3-4GBs mostly because I am exporting data as .MAT files)

What suggestions could be provided to improve the speed for such a stiff equation?

I have also tried stopping tests and it doesn't quite help all the time as I'd rather mathematica quit my program naturally as a result of overbearing stiffness than artificially through a stopping test. (The former has physical significance)

Yes, this question bears resemblance to this but I don't think its the same.

I have given Parallelize a thought but it doesn't work on NDSolve. Any options that I have either on the Mathematica front, computing front, or saving the interpolation function data?

# Some observation with LaunchKernel

## Edit:

Using the LaunchKernel[n] option just before the NDSolve cell didn't do much. My AbsoluteTiming barely even changed.

CloseKernels[];
LaunchKernels;
L = 2*92.389; TMax = 3100*100;
.........
......


## Edit 2:

By using Table and launching up to 6 kernels, these are the results that I got:

{{1,{19.454883,InterpolatingFunction[{{0.,184.778},{0.,184.778},{0.,282761.}},<>]}}, {2,{19.162008,InterpolatingFunction[{{0.,184.778},{0.,184.778},{0.,282761.}},<>]}}, {3,{18.919101,InterpolatingFunction[{{0.,184.778},{0.,184.778},{0.,282761.}},<>]}}, {4,{20.166785,InterpolatingFunction[{{0.,184.778},{0.,184.778},{0.,282761.}},<>]}}, {5,{20.265163,InterpolatingFunction[{{0.,184.778},{0.,184.778},{0.,282761.}},<>]}}, {6,{20.556365,InterpolatingFunction[{{0.,184.778},{0.,184.778},{0.,282761.}},<>]}}}

So with more kernels, the performance actually degraded....? Wha...?

• @NasserM.Abbasi Thanks for your comment. 1) How do you use the Interpolating function polynomial as an initial condition? 2) Define small time step – dearN Oct 6 '12 at 3:34
• @drN A "small time step" is something short like this, see? – Dr. belisarius Oct 6 '12 at 3:52
• @belisarius Did,you forget to link a page in your comment by any chance? :P` – dearN Oct 6 '12 at 13:46
• @NasserM.Abbasi Thanks! Will try that. I have a new issue though, is there any work around for "large" cases with lots of grid points, needed extensive amounts of RAM (8GB+)? – dearN Oct 6 '12 at 13:46