# 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\) +
Div[-h^3 Bo Grad[h] +
h^3 Grad[Laplacian[h]] + (\[Delta] h^3)/(Bi h + K1)^3 Grad[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
][[1]] // 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[3];
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 –  drN Oct 6 '12 at 3:34
@drN A "small time step" is something short like this, see? –  belisarius Oct 6 '12 at 3:52
@belisarius Did,you forget to link a page in your comment by any chance? :P` –  drN 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+)? –  drN Oct 6 '12 at 13:46
add comment

## 1 Answer

Yes, it is stiff -- but the main issue that I see is that the solution goes wild near the TMax that you specify. That's because you need a super-fine spatial grid to accurately represent what happens when the higher order terms finally manifest themselves. It's going to take a lot of time and a lot of memory (MinPoints option), and there's no way around it.

-
what about Remote Kernels? Any idea? I've posted an other question here –  drN Oct 7 '12 at 3:57
Nope. NDSolve isn't parallelized with Remote Kernels, which is why you got essentially the same 6 solutions in your previous post. However, it is parallelized with OpenMP for certain operations in certain solvers. There must be some linear algebra, FFTs, or matrix operations that they have multithreaded. You can see it if you monitor the CPU usage. –  Eric Brown Oct 7 '12 at 4:16
So, I am out of luck, eh? –  drN Oct 7 '12 at 4:23
I don't know, depends on how much time and memory you have to crank away with a larger grid. But if the calculation is mission critical, then you may want to try. Your problem is at least three-fold: 1) you probably need a specialized PDE solver (think Fortran, MPI, and big iron) due to the size of the grid, 2) it will still take a lot of time and memory, and 3) I would not like to reduce that 4th-order PDE down to first order to use in one of them. –  Eric Brown Oct 7 '12 at 4:42
Another thing to add is that you might want to try imtek but I have no experience with it. (The author's reputation is a first-rate Mathematica guru.) –  Eric Brown Oct 7 '12 at 4:46
add comment