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

Backslide introduced in 9.0, persisting through 11.3.

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.

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

• Seen this? – J. M. is away Feb 22 '16 at 5:46
• Thank you! This seems like what I am looking for! According to Prof. Tao, earthquake should be the key and therefore special initial or boundary conditions needed for the tsunami to propagate at the shore. I hope such may have been also discussed. Thank you. I am now reading it and try the notebook... – LCFactorization Feb 22 '16 at 5:51
• Related, but no nonlinearity: Modelling the effect of a structure on a “tsunami” (hyperbolic wave equation) – Jens Feb 22 '16 at 6:40
• About the solution of that linear partial differential equation, there is an example in Wolfram's NDSolve tutorial document here: reference.wolfram.com/language/tutorial/… , but is rather different from that introduced by Prof. Terence Tao – LCFactorization Feb 22 '16 at 10:35
• Can you show your previous efforts and add data/details ? At the moment, this sounds like a pure coding request... – Yves Klett Feb 24 '16 at 6:01

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

Go back to v8. v7 and v6 probably works, too.

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. I fail to find a simple workaround.

Another sad news is, if one increases nxy to e.g. 100, even v8 won't finish calculating in 500 seconds. We need to turn to an advanced approach to obtain a resulting graphic with higher resolution.

This problem turns out to be the 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 the 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 = DeveloperToPackedArray@
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 = CompileGetElement},
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} = NDSolveProcessEquations[{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"}}];

derifunc = Cases[state,a : NDSolveFiniteDifferenceDerivativeFunction[__][_] :> 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?

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 :) .