The Mathematics of Tsunamis

Mathematica code by Michael Trott
January, 2005
see http://mathworld.wolfram.com/news/2005-01-14/tsunamis/

Solve the radial shallow water equations using NDSolve

Here, we solve the shallow water wave equations for the radial case assuming an initial Gaussian water elevation in an “ocean” of constant depth in the middle and a cosine-like slope toward the shore.

((∂u(x,y,t))/(∂t))+u(x,y,t)((∂u(x,y,t))/(∂x)) +v(x,y,t)((∂u(x,y,t))/(∂y)) +g ((∂h(x,y,t))/(∂x))==0
((∂v(x,y,t))/(∂t))+u(x,y,t)((∂v(x,y,t))/(∂x)) +v(x,y,t)((∂v(x,y,t))/(∂y)) +g ((∂h(x,y,t))/(∂y))==0
((∂h(x,y,t))/(∂t))+((∂(u(x,y,t) (h(x,y,t)-b(x,y,t))))/(∂x))+((∂(v(x,y,t) (h(x,y,t)-b(x,y,t))))/(∂y))==0

Turn off spelling messages.

Off[General::"spell1"]

Set up the velocity and height functions.

u[x_, y_, t_] := Cos[φ] V[r, t]
v[x_, y_, t_] := Sin[φ] V[r, t]
h[x_, y_, t_] := h[r, t]
b[x_, y_, t_] := b[r]

Set up the shallow water wave partial differential equations.

eqs = 
Hold[
 {D[u[x, y, t], t] + u[x, y, t] D[u[x, y, t], x] + v[x, y, t] D[u[x, y, t], y] +
  g D[h[x, y, t], x] == 0,
  D[v[x, y, t], t] + u[x, y, t] D[v[x, y, t], x] + v[x, y, t] D[v[x, y, t], y] +
  g D[h[x, y, t], y] == 0, 
  D[h[x, y, t], t] + D[u[x, y, t] (h[x, y, t] - b[x, y, t]), x] +
                     D[v[x, y, t] (h[x, y, t] - b[x, y, t]), y] == 0}] /.
     {HoldPattern[D[a_, x]] :> Cos[φ] D[a, r], 
      HoldPattern[D[a_, y]] :> Sin[φ] D[a, r]};

Simplify them.

eqs1 = ReleaseHold[eqs] // FullSimplify

Use the simplified equations.

eqs = 
{(Derivative[0, 1][V][r, t] + g*Derivative[1, 0][h][r, t] + 
     V[r, t]*Derivative[1, 0][V][r, t]) == 0, 
  Derivative[0, 1][h][r, t] + V[r, t]*(-Derivative[1][b][r] + 
      Derivative[1, 0][h][r, t]) + (-b[r] + h[r, t])*
     Derivative[1, 0][V][r, t] == 0}

Define the bottom geometry.

b[r_] := -10 Piecewise[{{1, Sqrt[r^2] < R/3}, 
                        {Cos[1/2 Pi 3/2 (Sqrt[r^2] - R/3)/R]^2, Sqrt[r^2] < R}}]

Solve the equations using NDSolve (the real work).

ppR = 201;
g = 10; R = 10; RM = 5/4 R; T = 1;
ndsol = NDSolve[Join[eqs,
             {V[r, 0] == 0, h[r, 0] == Exp[-r^2],
              V[-RM, t] == 0, V[+RM, t] == 0}],
         {V, h}, {r, -RM, RM}, {t, 0, T},  
         MaxSteps -> 10^3, PrecisionGoal -> 3, AccuracyGoal -> 3, 
         Method -> {"MethodOfLines", 
                    "Method" -> "StiffnessSwitching",
                    "SpatialDiscretization" -> 
                    {"TensorProductGrid", "DifferenceOrder" -> Pseudospectral, 
                     "MaxPoints" -> {ppR}, "MinPoints" -> {ppR}}}]

Define some visualization routines.

rotationSurface[Line[l_]] := 
Module[{ℛ, points, ppφ = 60},
       points = Table[ℛ = {{Cos[-φ], Sin[-φ], 0}, {-Sin[-φ], Cos[-φ], 0},
                           {0, 0, 1}} // N;
                      ℛ.#& /@ l, {φ, 0, Pi, Pi/ppφ}]; 
          Table[Polygon[{points[[i, j]], points[[i + 1, j]], 
                         points[[i + 1, j + 1]], points[[i, j + 1]]}],
             {i, Length[points] - 1}, {j, Length[points[[1]]] - 1}]]

model3D[t_,opts___] := 
Module[{pp = 120, bottom, water, front, lB, lW},
 bottom = rotationSurface[Line[lB = N @ Table[{r, 0, b[r]}, {r, 0, R, R/pp}]]];
 water = rotationSurface[Line[lW = N @ Table[
                 Evaluate[{r, 0, h[r, t]} /. ndsol[[1]]], {r, 0, R, R/pp}]]];
 front = {#, Map[{-1, 1, 1}#&, #, {-2}]}&[Polygon[Join[#1, Reverse[#2]]]& @@@ 
                  Partition[Transpose[{lB, lW}], 2, 1]];               
Show[Graphics3D[{EdgeForm[], 
{{SurfaceColor[RGBColor[0, 1, 0], Hue[0.12], 2.2], bottom},
 {SurfaceColor[GrayLevel[0.7], RGBColor[0, 1, 0], 2.2], water, front}}}],
 opts, BoxRatios -> {1, 1, 0.3}, Axes -> False, Boxed -> False,
 ViewPoint -> {2, -3, 0.6}, PlotRange->{All,All,{-10.1,1.1}}]]

Make the animation.  (Note: this can take a while to render.)

tMax = 0.78; frames = 36;
movie = Table[model3D[t,ImageSize->500], {t, 0, tMax, tMax/frames}]