Introduction
A lack of time to write up an answer ironically provided time to reflect on the problem, and some nagging uncertainties about some issues contributed to the delay.
The slowness of the OP's code can be seen by a simple analysis, which reveals that some expensive calculations are repeated multiple times for each step in the time integration. A more important concern is whether area is conserved by the DE. The OP asserts yes, but I'm thinking not. I wonder which of us is correct. This issue has slowed my response, because I assumed area was conserved, which led to things not making sense. I first thought I could explain the change in area in terms of the finite-difference discretization. It does contribute to it, but not as much as the underlying differential equation. This was confusing to me until I accepted the hypothesis that area is not conserved. Similarly, a difference of opinion about the proper value of the grid spacing ds has made me hesitate some.
Code analysis
The OP's code for dh[] is clearly written and this is plugged into a simple straight-forward NDSolve call, so the whole is not hard to analyze.
Dependence and independence
I have found that I save myself some headaches if variables that depend on one another are constructed in code that reflects the dependency. Let's look at some of the quantities in the code, in particular s and ds.
First, ds and s are intimately linked to each, but their codes are independent. This is a weakness and a potential source for bugs. We should have ds = s[[2]] - s[[1]] (since it is a regular grid).
Similarly, h0 depends on n but is not defined in terms of n; indeed, it probably ought to be a function of θ and applied to s.
Control flow
The principal thing to notice is that in dh[], each of d1[#] and d2[#] is repeated several times and each entails a ListCorrelate[] computation to be repeated. This will slow things down unnecessarily.
By contrast, @xzczd's fast solution calculates dh[h] just once; even though it is a symbolic computation and somewhat slower than a numerical one, it is not nearly as slow as repeating calculations again and again throughout the steps of the time integration.
Leveraging Mathematica
1.
As shown in the tutorial, The Numerical Method of Lines, NDSolve`FiniteDifferenceDerivative[] will compute finite difference derivatives. It is virtually the same as the OP's ListCorrelation[] codes (with the appropriate "DifferenceOrder", with one difference: FiniteDifferenceDerivative[] uses a higher order approximation at the end points of the grid instead of padding (as the OP does with {2.}). The OP's results may be obtained by manually padding in the same way and dropping the first and last results from FiniteDifferenceDerivative[]. For instance, these give the same results as d1 and d2 respectively:
s = Array[# &, n, {-5/6 π, -π/6}];
ds = s[[2]] - s[[1]]; (* == Pi/300 for n = 201 *)
paddedgrid = Join[{First[s] - ds}, s, {Last[s] + ds}];
paddedh0 = Join[{2.}, h0, {2.}];
NDSolve`FiniteDifferenceDerivative[1, paddedgrid, (* == d1[h0] *)
"DifferenceOrder" -> 2][paddedh0][[2 ;; -2]]
NDSolve`FiniteDifferenceDerivative[2, paddedgrid, (* == d2[h0] *)
"DifferenceOrder" -> 1][paddedh0][[2 ;; -2]]
Note that NDSolve`FiniteDifferenceDerivative does not depend on ds and can be used just in terms of s.
2.
A minor point: The following is equivalent to the OP's use of ListCorrelation, padding h with 2. on the left and right:
d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, h, {-2, 2}, 2.]
d2 = ListCorrelate[{1, -2, 1}/ds^2, h, {-2, 2}, 2.]
3.
In addition, the OP's use of dh[h[t]], with dh[h_List], means that NDSolve will treat dh[] as a numeric black box, whereas , @xzczd's dh[h] constructs a (large) system of symbolic equations. It seems clear, from the difference in the number of steps shown below, that this affects the integration method used by NDSolve.
Improvement of speed in dh[]
Aside from eliminating the redundant computations outline above, the use of machine numbers, packed arrays, and vectorized operations, as discussed in
Performance tuning in Mathematica?, can be used to increase speed. Compile also helps quite a bit.
(The code for Compile was constructed from the OP's ode, and I used $ variables to avoid conflicting with other things.) One could do it with the padded grid shown above, which would add a little to the computation time, but I'll do it with the simpler approach.
ClearAll[dh2, s];
n = 201;
s = Developer`ToPackedArray@N@Array[# &, n, {-5/6 π, -π/6}];
dhC = Compile[{{$h, _Real, 1}, {$d1, _Real, 1}, {$d2, _Real, 1},
{$s, _Real, 1}, {$c, _Real, 1}},
-2 $d1 (-2 + $h)^2 (-1 + $h) ($c $h + $d1 $s) -
2 $d1 (-2 + $h) (-1 + $h)^2 ($c $h + $d1 $s) -
(-2 + $h)^2 (-1 + $h)^2 (2 $c $d1 + $d2 $s - $h $s)
];
With[{ (* compute derivatives and Sin[s], Cos[s] just once *)
d1FN = NDSolve`FiniteDifferenceDerivative[1, s, "DifferenceOrder" -> 2],
d2FN = NDSolve`FiniteDifferenceDerivative[2, s, "DifferenceOrder" -> 1],
sin = Sin[s],
cos = Cos[s]},
dh2[h_List] := With[{d1 = d1FN@h, d2 = d2FN@h},
dhC[h, d1, d2, sin, cos]]
];
Performance comparison
I wanted to compare evaluations of the right-hand side of the ODE, steps, and time. I added an option EvaluationMonitor :> ++evals to the NDSolve code for each case. I timed the solutions with code like the following one for my case:
{sol2} = h /. NDSolve[{h'[t] == dh2[h[t]], h[0] == h0}, h, {t, 0, 10},
EvaluationMonitor :> ++evals]; // AbsoluteTiming
Oddly, adding EvaluationMonitor to @xzczd's method increased the time by a huge factor (it took twice as long as mine), so I retimed it without the monitor, where @xzczd's took half the time mine did. These were done with a corrected ds.
One can see that dh2 evaluates about twice as fast, counted per evaluation, than @xzczd's method and 100 times faster than the OP's. Whether that is an entirely fair comparison depends on how much evaluating the ODE accounts for the execution time of NDSolve. It should account for a lot of it, but that's just my guess.
Area preserving?
I figure I should share this. If I made a mistake, well, someone can point it out. If not, then, maybe it will help the OP. If we consider what it means to be area-preserving, we should consider the integral
$$A = \int_{-\pi}^\pi \text{$1\over2$}\, h(t,\theta)^2\;d\theta \,,$$
where $h$ is extended to $-\pi \le \theta \le \pi$ by defining $h(t,\theta) = 2$ (that is, the greater radius) outside $-5\pi/6 \le \theta \le -\pi/6$.
Then, by my lights, differentiating under the integral sign, the rate of change of the area is given by
$${dA \over dt} = \int_{-\pi}^\pi h(t,\theta)h_t(t,\theta)\;d\theta$$
and $h_t$ can be written in terms of $h$, $h_\theta$, $h_{\theta\theta}$ using the PDE.
Of course, an equilibrium solution will preserve area. But if we perturb an equilibrium solution, then $dA/dt$ is nonzero:
ClearAll[dhdt, dAdt];
Block[{h, θ},
dhdt[h_, θ_] = -D[(h[θ] - 2)^2 (h[θ] - 1)^2 D[h[θ] Sin[θ], θ], θ];
];
dAdt[hfn_] := NIntegrate[hfn[θ] dhdt[hfn, θ], {θ, -Pi, Pi}];
(* background *)
θ1 = -2.4539664700963293`;
background = Show[
RegionPlot[1 <= x^2 + y^2 <= 4, {x, -2, 2}, {y, -2, 2},
PlotStyle -> Opacity[0.1], AspectRatio -> 1],
PolarPlot[
2 Sin[θ1]/Sin[θ], {θ, θ1, -(Pi + θ1)},
PlotStyle -> Black]
];
solPlot[hfn_] := Show[
background,
PolarPlot[hfn, {θ, -5/6 π, -π/6},
ColorFunction -> (ColorData["RedGreenSplit"][
0.5 (1 + Clip[100 dhL /. h -> Function @@ {θ, hfn} /. θ -> #3, {-1,1}])] &),
ColorFunctionScaling -> False,
PlotRange -> {{-2, 2}, {-2, 2}}]
];
Block[{Print, θ2 = θ1, eps = +0.1, hf},
hf = Piecewise[{
{2 Sin[θ2] Csc[#] + eps (# - θ2) (# + Pi + θ2), θ2 < # < -Pi - θ2}
}, 2] &;
Show[solPlot[hf[θ]], PlotLabel -> Row[{"dA/dt = ", dAdt[hf]}]]
]
(Red means the radius h is decreasing; green means it is increasing (with respect to time in the DE.)
The effect of finite differences and the error in ds
Using Block and the OP's definition of dh, we can see what happens when we use different values for ds. There is some controversy about the correct value for ds, and it is my view that it should be equal to the differences in s. The padding by 2 on either side of h in for ListCorrelate[] in the OP's dh[] corresponds to using s == -5 Pi/ 6 - ds and s == -Pi / 6 + ds and corresponds to the (unused) condition that h is identically equal to 2 outside the interval {-5 Pi/6, -Pi/6}.
If we plot dh[h[s]], where h is the equilibrium solution of the OP, h[s] = 2 Sin[θ1] Csc[s] for a given θ1, and s is the input grid as above, we will see how close to zero the value of dh is for the finite difference scheme.
ClearAll[dhEquilPlot];
SetAttributes[dhEquilPlot, Listable];
θ1 = -2.4539664700963293`;
dhEquilPlot[ds0_] := Block[{ds = ds0},
ListPlot[
dh[Clip[2 Sin[θ1] Csc[s], {0., 2.}]],
PlotLabel -> Row[{"ds = ", ds}], Frame -> True,
FrameLabel -> {HoldForm[s], HoldForm[h'[t]]}, DataRange -> MinMax[s],
BaseStyle -> {FractionBoxOptions -> {Beveled -> True}}
]
]
dhEquilPlot /@ {
{1/n}, (* OP (UL) *)
{(-(π/6) - -((5 π)/6))/(n + 2 - 1), (* xzczd (LL) *)
(-(π/6) - -((5 π)/6))/(n - 1)} (* M E2 (LR) *)
} // GraphicsGrid
We can see how much closer to zero dh is for ds = Pi/300, which is also the grid spacing in s.
It is also clear that in all three cases, the radius h is pulled toward the center (dh < 0) near where the level of the fluid contacts the outside circle (where h reaches 2). Of course the actual evolution from the OP's starting configuration will be slightly different. But since physically, the model should be area-preserving, there is little point in exploring this beyond determining the value of ds until the issue with area is cleared up.
valu? The argument ofdh[h]is limited toList, sodh[h[t]]will just return unevaluated. – xzczd Jun 26 '16 at 09:30valu, it's meant to be avalh. Don't know what you mean aboutdh[h[t]]being unevaluated. TheNDSolvereturns anInterpolatingFunctionof dimension n as required? – mch56 Jun 26 '16 at 11:45NDSolvefor handling vectors is improved, but didn't expect that it can handle function definition likeh[t_List]now. Just a side note: in a 2GHz laptop it takes about 85 seconds to solve for{t, 0, 1}. – xzczd Jun 26 '16 at 12:42h[-Pi/6] ==h[-5Pi/6]==2? – xzczd Jun 27 '16 at 10:15tto something like{t, 0, 20}(This takes about 140 seconds in my PC, of course you needMaxSteps->Infinityoption in this case) and observe the result, you'll see the solution does converge to a steady state, but not the one you're expecting. – xzczd Jun 27 '16 at 11:16NDSolvewithout writing my own FD scheme. I observe that usingNDSolvethe PDE tends to the correct (physically intuitive) steady state. If I then use this same initial condition in my FD scheme, the PDE tends to the wrong steady state. Therefore I conclude there must be an error with my scheme. It's a shameNDSolvecan't handle my desired initial condition without getting itself in a twist. I can provide further details on the working smooth initial condition case if you'd like? – mch56 Jun 27 '16 at 11:37NDSolvecan't handle my desired initial condition without getting itself in a twist." So a fullyNDSolvebased solution is OK to you, ifNDSolveis able to solve it? As far as I can tell,NDSolvecan handle discontinued initial condition (at least in some cases), how did you set that discontinued initial condition inNDSolve? – xzczd Jun 27 '16 at 11:53h0.NDSolvegets upset about the equation being stiff. This is because of the shock formation near the nose of the current - you can see this in the first gif. I would happily useNDSolveif it could handle these large gradients, but I'm not sure how? This is why I decided to write my own FD instead. – mch56 Jun 27 '16 at 13:02