Skip to main content

differential equations - 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["DifferentialEquations`InterpolatingFunctionAnatomy`"];
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?



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...?



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.


Comments

Post a Comment

Popular posts from this blog

plotting - Plot 4D data with color as 4th dimension

I have a list of 4D data (x position, y position, amplitude, wavelength). I want to plot x, y, and amplitude on a 3D plot and have the color of the points correspond to the wavelength. I have seen many examples using functions to define color but my wavelength cannot be expressed by an analytic function. Is there a simple way to do this? Answer Here a another possible way to visualize 4D data: data = Flatten[Table[{x, y, x^2 + y^2, Sin[x - y]}, {x, -Pi, Pi,Pi/10}, {y,-Pi,Pi, Pi/10}], 1]; You can use the function Point along with VertexColors . Now the points are places using the first three elements and the color is determined by the fourth. In this case I used Hue, but you can use whatever you prefer. Graphics3D[ Point[data[[All, 1 ;; 3]], VertexColors -> Hue /@ data[[All, 4]]], Axes -> True, BoxRatios -> {1, 1, 1/GoldenRatio}]
Post a Comment
Read more

plotting - Filling between two spheres in SphericalPlot3D

Manipulate[ SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, Mesh -> None, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], {n, 0, 1}] I cant' seem to be able to make a filling between two spheres. I've already tried the obvious Filling -> {1 -> {2}} but Mathematica doesn't seem to like that option. Is there any easy way around this or ... Answer There is no built-in filling in SphericalPlot3D . One option is to use ParametricPlot3D to draw the surfaces between the two shells: Manipulate[ Show[SphericalPlot3D[{1, 2 - n}, {θ, 0, Pi}, {ϕ, 0, 1.5 Pi}, PlotPoints -> 15, PlotRange -> {-2.2, 2.2}], ParametricPlot3D[{ r {Sin[t] Cos[1.5 Pi], Sin[t] Sin[1.5 Pi], Cos[t]}, r {Sin[t] Cos[0 Pi], Sin[t] Sin[0 Pi], Cos[t]}}, {r, 1, 2 - n}, {t, 0, Pi}, PlotStyle -> Yellow, Mesh -> {2, 15}]], {n, 0, 1}]
Post a Comment
Read more

plotting - Mathematica: 3D plot based on combined 2D graphs

I have several sigmoidal fits to 3 different datasets, with mean fit predictions plus the 95% confidence limits (not symmetrical around the mean) and the actual data. I would now like to show these different 2D plots projected in 3D as in but then using proper perspective. In the link here they give some solutions to combine the plots using isometric perspective, but I would like to use proper 3 point perspective. Any thoughts? Also any way to show the mean points per time point for each series plus or minus the standard error on the mean would be cool too, either using points+vertical bars, or using spheres plus tubes. Below are some test data and the fit function I am using. Note that I am working on a logit(proportion) scale and that the final vertical scale is Log10(percentage). (* some test data *) data = Table[Null, {i, 4}]; data[[1]] = {{1, -5.8}, {2, -5.4}, {3, -0.8}, {4, -0.2}, {5, 4.6}, {1, -6.4}, {2, -5.6}, {3, -0.7}, {4, 0.04}, {5, 1.0}, {1, -6.8}, {2, -4.7}, {3, -1....
Post a Comment
Read more