How to solve the tsunami model and animate the shallow water wave?

Question

Backslide introduced in 9.0, persisting through 13.1.

---

Recently when I was learning differential equations, I noticed there is a shallow water wave equation to model the tsunami propagation.

How to establish and solve the initial and boundary conditions for the partial differential equations system of the shallow water wave model by NDSolve? and how to animate the tsunami at shore using the model in Mathematica for demonstration purpose of the model?

Update

thanks to @J.M. 's comment, I find there is already a sample in a Wolfram blog by Rob Knapp.

The visualization by Rob Knapp in a video is very vivid. Now I would like to update the question to below so that any others can post an answer:

1. solve the current problem with variable initial and boundary conditions, or the current ibcs a little faster or not so demanding in memory;

It seems visualization of the ocean floor together with the wave thus obtained by Mathematica remains a challenge.

2. visualize the video similar to that appears in Rob Knapp's blog by Mathematica

The first answer that solves any of the two will be accepted.

Answer 1

The most hard part of OP's questions is the latter half of 1st one i.e.

solve the problem in a little faster or not so demanding-in-memory way.

I've found 2 solutions, one is easy but hard to extend, the other is advanced but general. Let me talk about the simple one first.

Simple solution

Manually set Method -> "Adams" or Method -> "LSODA" for time integration i.e. changing the option in NDSolve to

Method -> {"MethodOfLines", Method -> "Adams" (* Alternatively: "LSODA" *), 
           "SpatialDiscretization" -> 
              {"TensorProductGrid", "MaxPoints" -> nxy, "MinPoints" -> nxy, 
                "DifferenceOrder" -> "Pseudospectral"}}

will speed up the calculation in or before v8, and seems to be necessary to make NDSolve solve out the system since v9. Clearly there exists kind of backslide between v8 and v9. I've tested the code in v9.0.1 and v11.2, but none of them manages to finish calculating in 500 seconds without Method -> "Adams" or Method -> "LSODA".

Remark

Method tracing with the approach in this post and time step monitoring with the approach in this post suggests NDSolve uses "StiffnessSwitching" method to solve the system by default, but still, what change has happened between v8 and v9 remains unclear.

With Method -> "Adams" or Method -> "LSODA", on a 2GHz dual core machine, the timing of NDSolve is about 11 seconds for nxy = 50, and about 80 seconds for nxy = 100. Not too bad, but can be faster if we turn to the following advanced approach.

Advanced solution

This problem turns out to be 3rd one I find in this site, on which the default spatial discretization of NDSolve doesn't work well even if nothing like shock wave involves in. (BTW here are 1st one and 2nd one. ) If you set e.g. "DifferenceOrder" -> 2 instead of "DifferenceOrder" -> "Pseudospectral", you'll find the solution becomes unstable soon. The following is the solution for nxy = 50, "DifferenceOrder" -> 2 at t == T:

However, 2nd order difference scheme is enough to solve this problem, if only derivative terms like $\frac{\partial (f g)}{\partial x}$ is discretized as

$$\frac{\partial (f g)}{\partial x}\Biggl|_{x=x_i} \approx \frac{f(x_i+h) g(x_i+h)-f(x_i-h)g(x_i-h)}{2h}$$

Notice this scheme is different from:

$$\frac{\partial (f g)}{\partial x}\Biggl|_{x=x_i}=\left(f\frac{\partial g}{\partial x}+g\frac{\partial f}{\partial x}\right)\Biggl|_{x=x_i}\approx f(x_i) \frac{g(x_i+h)-g(x_i-h)}{2h}+g(x_i)\frac{f(x_i+h)-f(x_i-h)}{2h}$$

which is the one used by NDSolve when "DifferenceOrder" -> 2.

Remark

Though "DifferenceOrder" and actual difference order chosen by NDSolve happen to be the same in this case, their relationship isn't that simple. Check this post for more information.

The remaining work is coding:

Seamount[s_, w_, {λ_, ϕ_}][l_, p_] := s Exp[-w ((l - λ)^2 + (p - ϕ)^2)]; 
FaultDisplacement[p1_, p2_][p_?(VectorQ[#, NumberQ] &)] :=
 If[p1 == p2, Norm[p - p1],
   Module[{pp1, pp2, p12},
   pp1 = p - p1;
   pp2 = p - p2;
   p12 = (p1 - p2)/Norm[p1 - p2];
   If[Sign[p12.pp1] != Sign[p12.pp2],
    Norm[pp1 - pp1.p12 p12],
    Min[Norm[pp1], Norm[pp2]]]]];

λmin = -15 ° ; λmax = 15 °; λe = 0; λi = 2.5 °;
ϕmin = -35 °; ϕmax = 5 °; ϕ1 = -10 °; ϕ2 = -20 °; ϕi = -25 °;
b[λ_, ϕ_] := -4000 + 
   Seamount[3500, 1000, {-2.5 °, -25 °}][λ, ϕ] + 
   Seamount[3000, 1000, {-5 °, -23 °}][λ, ϕ] +  
   Seamount[3000, 1000, {-7.5 °, -27 °}][λ, ϕ] +  
   Seamount[3000, 1000, {5 °, -24 °}][λ, ϕ] + 
   Seamount[3000, 1000, {5 °, -21 °}][λ, ϕ] + 
   Seamount[3000, 1000, {5 °, -18 °}][λ, ϕ] + 
   Seamount[3000, 1000, {5 °, -15 °}][λ, ϕ] + 
   Seamount[3000, 1000, {5 °, -12 °}][λ, ϕ] +
   Seamount[3750, 1000, {-2 °, -18 °}][λ, ϕ];
ClearAll[ct]
SetAttributes[#, HoldAll] & /@ {ct};
ct@D[expr_, x_] :=Subtract@@(expr /. {{x -> x +delta@x}, {x -> x - delta@x}})/(2 delta@x)

points@λ = 300; points@ϕ = 300;
delta@λ = (λmax - λmin)/(points@λ - 1);
delta@ϕ = (ϕmax - ϕmin)/(points@ϕ - 1);
With[{h = h[λ, ϕ], u = u[λ, ϕ], b = b[λ, ϕ], v = v[λ, ϕ]}, 
  rhs= Unevaluated@
      {-(1/(R Cos[ϕ])) (D[u (h - b), λ] + D[v (h - b) Cos[ϕ], ϕ]), 
       -((D[u, ϕ]  v )/R + (g D[h, λ])/(R Cos[ϕ]) + (u  D[u, λ])/(R Cos[ϕ])) + f v,
      -((D[v, ϕ]  v )/R + (g D[h, ϕ])/R + (u  D[v, λ])/(R Cos[ϕ])) - f u 
      } /. D -> (ct@D[##] &)];
iclst = Developer`ToPackedArray@
   Transpose[Table[
     N@{ Exp[-2000 FaultDisplacement[{λe, ϕ1}, {λe, ϕ2}][{λ, ϕ}]^2], 0, 0}, 
      {λ, λmin, λmax, delta@λ}, {ϕ, ϕmin, ϕmax, delta@ϕ}], {2, 3, 1}];
rt = RescalingTransform[{{λmin, λmax}, {ϕmin, ϕmax}}, {{1, points@λ}, {1, points@ϕ}}];

{R = 6378000, ω = (2 π)/(24 3600), f = 2 ω Sin[ϕ], g = 9.8, s = .95};
With[{rc = RuleCondition, cg = Compile`GetElement}, 
  rhsfunc = Hold@
           Compile[{{var, _Real, 3}}, expr, 
            RuntimeOptions -> EvaluateSymbolically -> False, CompilationTarget -> C] /. 
          expr -> (table[#, {λ, λmin, λmax, delta@λ}, {ϕ, ϕmin, ϕmax, delta@ϕ}] & /@ 
             rhs) /. table -> Table /. Thread[{h, u, v} -> (var /@ Range@3)] /. 
       var[i_][lambda_, phi_] :> 
        rc@(cg[var, i, Mod[#, points@λ - 1, 1], 
             Mod[#2, points@ϕ - 1, 1]] & @@ Round /@ rt@{lambda, phi}) /. 
      Degree -> N@Degree // ReleaseHold // Last];
T = 2 3600; 
 resultfunclst = 
  NDSolve[{var'[t] == rhsfunc@var@t, var[0] == iclst}, var, {t, 0, T}, 
    MaxSteps -> Infinity][[1, 1, -1]]; // AbsoluteTiming
(* {334.9169883,Null} *)

Remark

It's not immediately clear what boundary conditions (b.c.s) are chosen for u and v in the original code, but after some spelunking I found the b.c.s for u and v are also periodic b.c.:

Block[{R = 6378000, ω = (2 π)/(24 3600), f = 2 ω Sin[ϕ], g = 9.8, s = .95},
 With[{h = h[λ, ϕ, t], u = u[λ, ϕ, t], b = b[λ, ϕ], v = v[λ, ϕ, t]}, 
  swe = {D[h, t] + 1/(R Cos[ϕ]) (D[u (h - b), λ] + D[v (h - b) Cos[ϕ], ϕ]) == 0, 
         (D[u, ϕ] v)/R +(g D[h, λ])/(R Cos[ϕ]) + D[u, t] + (u D[u, λ])/(R Cos[ϕ])==f v,
         (D[v, ϕ] v)/R + (g D[h, ϕ])/R + D[v, t] + (u D[v, λ])/(R Cos[ϕ]) == -f u};
  ic = {h == Exp[-2000 FaultDisplacement[{λe, ϕ1}, {λe, ϕ2}][{λ, ϕ}]^2], 
        u == 0, v == 0} /. t -> 0;
  bc = {(h /. ϕ -> ϕmin) == (h /. ϕ -> ϕmax), 
        (h /. λ -> λmin) == (h /. λ -> λmax)}]];

nxy = 50;
{state} = NDSolve`ProcessEquations[{swe, ic , bc}, 
   h, {λ, λmin,  λmax}, {ϕ, ϕmin, ϕmax}, {t, 0, 
    T}, DependentVariables -> {h, u, v}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> nxy, 
       "MinPoints" -> nxy, "DifferenceOrder" -> "Pseudospectral"}}];

ml = LinkCreate[LinkMode -> Loopback];
LinkWrite[ml, With[{e = state}, Hold[e]]];
holdstate = LinkRead[ml];

derifunc = Cases[If[$VersionNumber > 11.2, holdstate, state], 
   a : NDSolve`FiniteDifferenceDerivativeFunction[__][_] :> Head@a, Infinity] // Union;

#["PeriodicInterpolation"] & /@ derifunc
(* {{True, True}, {True, True}, {True, True}, 
    {True, True}, {True, True}, {True, True}} *)

So periodic b.c. is set in both directions for all 3 variables in my code.

If you don't have a C compiler installed, take the CompilationTarget option and // Last away. I do recommend you to install one though.

Some advanced technique has been used in this code piece to make the discretization less tedious. To understand it better, you may want to read the following posts:

When should I, and when should I not, set the HoldAll attribute on a function I define?

How to make the code inside Compile conciser without hurting performance?

Why is CompilationTarget -> C slower than directly writing with C?

Replacement inside held expression

Notice I've chosen a dense grid here. With points@λ = 50; points@ϕ = 50; NDSolve will finish calculating in about 2 seconds and the max memory used is about 50 MB. (Of course the error will be a bit large, but not too bad in my view. )

With the advanced solution at hand, 2nd question is trivial:

oceanfloor = 
 Plot3D[b[λ, ϕ] // 
   Evaluate, {λ, λmin, λmax}, {ϕ, ϕmin, ϕmax}, 
  ColorFunction -> (ColorData["GreenBrownTerrain"][0.05 + 0.95 #3] &), PlotRange -> All, 
  Mesh -> None, PlotPoints -> 50]

sf = 3 10^3; 
lst = Table[
   ListPlot3D[resultfunclst[t][[1]]\[Transpose] sf, Mesh -> None, 
     DataRange -> {{λmin, λmax}, {ϕmin, ϕmax}}, 
     PlotStyle -> Opacity[1/2], PlotRange -> {-4000, 3000}, Axes -> None, 
     Boxed -> False]~Show~oceanfloor, {t, 0, T, 360}];

ListAnimate[lst]

MaxMemoryUsed[]/1024^2.
(* 2965.88 *)

Though I fail to reproduce the texture of sea floor, I think the animation is similar enough to the one in the movie.

"OK, then how about the first half of 1st question? " Well, modifying the initial condition (i.c.) in the advanced solution is also trivial, but making the code accept other type of b.c. is not that simple. (For b.c. involving derivative, technique like one-sided difference formula is needed. ) Given the PDE system is for simulating tsunami, I think the current b.c. i.e. periodic b.c. approximating b.c. at infinity is quite enough, so I'd like to skip this request. Anyway, I've already overfulfilled :) .