Shallow-water equations on a rotating sphere by using the NDSolve method with a purpose of earthquake-generated tsunami simulation

Question

I'm trying to numerically solve a system of Shallow Water Equations in Wolfram Mathematica 14.0.0 on the rotating sphere by using the NDSolve method with a purpose of earthquake-generated Tsunami simulation. I define necessary physical constants, initial ocean floor profile, initial values of angular speeds and water surface. I also setup time-dependent ocean floor excitation law and atmospheric pressure function.

Shallow Water Gravity Waves Equations are presented below: Here are my code cuttings:

r0 = 6.4*10^6; (*Earth radius*)
g = 9.81; (*Earth gravitational constant*)
\[Omega]0 = 7.3*10^-5; (*Earth angular speed*)
\[Omega]0 = 0; (*set in order to simplify problem visualization*)
d = 10^4; (*average depth of ocean between surface and floor*)
p0 = 1.013*10^5; (*normal atmospheric pressure*)
\[Rho] = 10^3; (*water normal density*)
\[Nu] = 1.003*10^-6; (*water normal kinematic viscosity*)

\[Theta]min = 0 + \[Pi]/10;
\[Theta]max = \[Pi] - \[Pi]/10;
\[Phi]min = 0;
\[Phi]max = 2 \[Pi];

FloorProfile[\[Theta]_, \[Phi]_] := 
  d/4*Sin[2 \[Theta]]*(1 + 1/2 Sin[\[Phi]]);
Compensator[\[Theta]_, \[Phi]_] := 
  1/(2*g)*(\[Omega]0*r0)^2*Sin[\[Theta]]^2;

\[Theta]exmin = 0.45*\[Pi];
\[Theta]exmax = 0.55*\[Pi];
\[Phi]exmin = 0.95*\[Pi];
\[Phi]exmax = 1.05*\[Pi];
aex = d/100;
tex = 10;

FloorProfileExcitation[t_, \[Theta]_, \[Phi]_] := 
  If[t <= 0, 0, If[t >= tex, aex, aex/tex*t]]*
   If[\[Theta]exmin < \[Theta] < \[Theta]exmax, 
    Sin[\[Pi]*((\[Theta] - \[Theta]exmin)/(\[Theta]exmax - \
\[Theta]exmin))]^2, 0]*
   If[\[Phi]exmin < \[Phi] < \[Phi]exmax, 
    Sin[\[Pi]*((\[Phi] - \[Phi]exmin)/(\[Phi]exmax - \
\[Phi]exmin))]^2, 0];

f[t_, \[Theta]_, \[Phi]_] := 
  FloorProfile[\[Theta], \[Phi]] + 
   FloorProfileExcitation[t, \[Theta], \[Phi]];
p[t_, \[Theta]_, \[Phi]_] := p0;
h0[\[Theta]_, \[Phi]_] := Compensator[\[Theta], \[Phi]];
v0[\[Theta]_, \[Phi]_] := 0;
u0[\[Theta]_, \[Phi]_] := 0;

tsim = 10^5;

equations = {r0*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\((h[t, \ \[Theta], \ \[Phi]] - 
         f[t, \ \[Theta], \ \[Phi]])\)\) + 
     v[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((h[
          t, \ \[Theta], \ \[Phi]] - f[t, \ \[Theta], \ \[Phi]])\)\) +
      u[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\((h[
          t, \ \[Theta], \ \[Phi]] - f[t, \ \[Theta], \ \[Phi]])\)\)/
      Sin[\[Theta]] + (h[t, \[Theta], \[Phi]] - 
        f[t, \[Theta], \[Phi]] + d)*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((v[
            t, \ \[Theta], \ \[Phi]]*Sin[\[Theta]])\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(u[
          t, \ \[Theta], \ \[Phi]]\)\))/Sin[\[Theta]] == 0, r0*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(v[
        t, \ \[Theta], \ \[Phi]]\)\) + v[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(v[
        t, \ \[Theta], \ \[Phi]]\)\) + u[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(v[
        t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] == 
    u[t, \[Theta], \[Phi]]^2*Cos[\[Theta]]/Sin[\[Theta]] + 
     2*\[Omega]0*r0*u[t, \[Theta], \[Phi]]*
      Cos[\[Theta]] + (\[Omega]0*r0)^2 Sin[\[Theta]]*Cos[\[Theta]] - 
     g*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(h[
        t, \ \[Theta], \ \[Phi]]\)\) - 1/\[Rho]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(p[
        t, \ \[Theta], \ \[Phi]]\)\) + \[Nu]/r0*( \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((Sin[\[Theta]]*\ 
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]v[
             t, \ \[Theta], \ \[Phi]])\)\)/Sin[\[Theta]] + \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]v[
           t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]]^2), r0*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(u[
        t, \ \[Theta], \ \[Phi]]\)\) + v[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(u[
        t, \ \[Theta], \ \[Phi]]\)\) + u[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(u[
        t, \ \[Theta], \ \[Phi]]\)\)/
      Sin[\[Theta]] == -v[t, \[Theta], \[Phi]]*u[t, \[Theta], \[Phi]]*
      Cos[\[Theta]]/Sin[\[Theta]] - 
     2*\[Omega]0*r0*v[t, \[Theta], \[Phi]]*Cos[\[Theta]] - g*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(h[
        t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] - 1/\[Rho]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(p[
        t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] + \[Nu]/r0*( \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((Sin[\[Theta]]*\ 
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]u[
             t, \ \[Theta], \ \[Phi]])\)\)/Sin[\[Theta]] + \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]u[
           t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]]^2)};


initialconditions1 = {h[0, \[Theta], \[Phi]] == h0[\[Theta], \[Phi]], 
   v[0, \[Theta], \[Phi]] == v0[\[Theta], \[Phi]], 
   u[0, \[Theta], \[Phi]] == u0[\[Theta], \[Phi]], 
   PeriodicBoundaryCondition[
    h[t, \[Theta], \[Phi]], \[Phi] == 2*\[Pi], 
    TranslationTransform[{0, -2*\[Pi]}]], 
   PeriodicBoundaryCondition[
    v[t, \[Theta], \[Phi]], \[Phi] == 2*\[Pi], 
    TranslationTransform[{0, -2*\[Pi]}]], 
   PeriodicBoundaryCondition[
    u[t, \[Theta], \[Phi]], \[Phi] == 2*\[Pi], 
    TranslationTransform[{0, -2*\[Pi]}]], 
   PeriodicBoundaryCondition[h[t, \[Theta], \[Phi]], \[Phi] == 0, 
    TranslationTransform[{0, 2*\[Pi]}]], 
   PeriodicBoundaryCondition[v[t, \[Theta], \[Phi]], \[Phi] == 0, 
    TranslationTransform[{0, 2*\[Pi]}]], 
   PeriodicBoundaryCondition[u[t, \[Theta], \[Phi]], \[Phi] == 0, 
    TranslationTransform[{0, 2*\[Pi]}]]}; 
method1 = 
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "TemporalVariable" -> t, 
      "SpatialDiscretization" -> {"FiniteElement", 
        "MeshOptions" -> {"MeshElementType" -> "TriangleElement"}}}};
region1 = {\[Theta], \[Phi]} \[Element] 
   Rectangle[{\[Theta]min, \[Phi]min}, {\[Theta]max, \[Phi]max}];

initialconditions2 = {h[0, \[Theta], \[Phi]] == h0[\[Theta], \[Phi]], 
   v[0, \[Theta], \[Phi]] == v0[\[Theta], \[Phi]], 
   u[0, \[Theta], \[Phi]] == u0[\[Theta], \[Phi]], 
   h[t, \[Theta], 0] == h[t, \[Theta], 2 \[Pi]], 
   v[t, \[Theta], 0] == v[t, \[Theta], 2 \[Pi]], 
   u[t, \[Theta], 0] == u[t, \[Theta], 2 \[Pi]]};
method2 = 
  Method -> {"PDEDiscretization" -> {"MethodOfLines", 
      "TemporalVariable" -> t, 
      "SpatialDiscretization" -> {"TensorProductGrid", 
        "DifferenceOrder" -> "Pseudospectral"}}};
region2 = 
  Sequence @@ {{\[Theta], \[Theta]min, \[Theta]max}, {\[Phi], \
\[Phi]min, \[Phi]max}};

initialconditions3 = {h[0, \[Theta], \[Phi]] == h0[\[Theta], \[Phi]], 
   v[0, \[Theta], \[Phi]] == v0[\[Theta], \[Phi]], 
   u[0, \[Theta], \[Phi]] == u0[\[Theta], \[Phi]]};
method3 = Method -> Automatic;

RunNDSolveValue[opts___] := 
  NDSolveValue[{equations, initialconditions1}, {h, v, u}, {t, 0, 
    tsim}, region1, MaxSteps -> 5000, 
   MaxStepSize -> 50, StartingStepSize -> 1, 
   InterpolationOrder -> All, AccuracyGoal -> Automatic, method1, 
   PrecisionGoal -> Automatic, opts];
Module[{time = 
     0}, (Monitor[
      solution = RunNDSolveValue[StepMonitor :> (time = t)], 
      ProgressIndicator[time/tsim]];)]; // AbsoluteTiming

My WM implementation code with plots and some other metadata was uploaded on my Google Drive.

In general, I need to simulate wave propagation on a very specific domain (world-wide ocean landscape) with periodic boundary conditions over latitude.

According to the post, my case could be evaluated only with FiniteElement method (FEM). As it was shown in @SantaP's answer under the post (thanks @xzczd for link), I ought to use PeriodicBoundaryCondition with "MeshOptions" -> {"MeshElementType" -> "TriangleElement"} option for FEM in order to obtain stable correct solution.

So, my first toy model is expressed via initialconditions1 & method1 & region1 parameters for RunNDSolveValue function.

However, at some point of simulation I always receive an error "At t == ~18000, step size is effectively zero; singularity or stiff system suspected." When I try to visualize my solution, I also see that it tends to rapidly burst after some iterations. One could see my problem on the animated video below:

Being frustrated about FEM's accuracy by the discussion, I've tried to play with "SpatialDiscretization" -> {"TensorProductGrid", "DifferenceOrder" -> "Pseudospectral"} as it was recommended in the post. This attempt corresponds to initialconditions2 & method2 & region2 parameters for RunNDSolveValue function. Boundary conditions are still periodic here, but further complexification of the spatial domain is unavailable. Unfortunately, this hopeful scheme also crashes with the same error "At t == ~1200, step size is effectively zero; singularity or stiff system suspected." and new warning "An insufficient number of boundary conditions have been specified for the direction of independent variable [Theta]. Artificial boundary effects may be present in the solution."

Moreover, the most trivial default technique initialized with initialconditions3 & method3 & region2 parameters for RunNDSolveValue has the similar behaviour.

Furthermore, I’ve also googled some other different options for NDSolveValue function from here and here, but none of them helped me to stabilize calculation process near the boundary of the simulation region. At the small simulation times I can clearly see my tsunami propagating correctly. In contrast, at the large times numerical interpolation solution explodes with enormous water surface height values without any reason.

So, give me, please, some advice how to deal with such issue. I’m out of ideas ☹

Answer 1

We can solve this problem using Method of Lines up to t=0.89 timeSim=89000. But realistic solution we can visualize up to t=30000. Next, we have problems with the wave reaching the gentle shore (North Pole). Note that I have changed the initial and boundary conditions. As a timescale we used timeSim.

 \[CapitalOmega]0 = 0 * 7.3*10^-5;
a0 = 6.4*10^6;
H0 = 10^4;
g = 10;
P0 = 101300;
\[Rho] = 1000;
minTheta = 0;
maxTheta = \[Pi];
minPhi = 0;
maxPhi = 2*\[Pi];
OceanProfile[\[Theta]_, \[Phi]_] := 
  H0/4*Sin[2 \[Theta]]*(1 + 1/2 Sin[\[Phi]]);
Compensator[\[Theta]_, \[Phi]_] := 
  1/(2*g)*(\[CapitalOmega]0*a0)^2*Sin[\[Theta]]^2;
minThetaEx = 0.45*\[Pi];
maxThetaEx = 0.55*\[Pi];
minPhiEx = 0.95*\[Pi];
maxPhiEx = 1.05*\[Pi];
timeSim = 10^5;
amplitudeEx = H0/100;
timeEx = 10/timeSim;
OceanProfileEx = 
  Compile[{{t, _Real}, {\[Theta], _Real}, { \[Phi], _Real}}, 
   Module[{}, 
    If[t >= timeEx, 0, amplitudeEx/timeEx *t]*
     If[minThetaEx < \[Theta] < maxThetaEx, 
      Sin[\[Pi]*((\[Theta] - minThetaEx)/(
         maxThetaEx - minThetaEx))]^2, 0]*
     If[minPhiEx < \[Phi] < maxPhiEx, 
      Sin[\[Pi]*((\[Phi] - minPhiEx)/(maxPhiEx - minPhiEx))]^2, 0]]];
profileFullFunction[t_, \[Theta]_, \[Phi]_] := 
  OceanProfile[\[Theta], \[Phi]] + OceanProfileEx[t, \[Theta], \[Phi]];
patmFullFunction[t_, \[Theta]_, \[Phi]_] := 0;
h0Function[\[Theta]_, \[Phi]_] := Compensator[\[Theta], \[Phi]];
v0Function[\[Theta]_, \[Phi]_] := 0;
w0Function[\[Theta]_, \[Phi]_] := 0;

Plot[OceanProfileEx[
  t, (maxThetaEx + minThetaEx)/2, (maxPhiEx + minPhiEx)/2], {t, 0, 
  10^-3}]

a = a0/timeSim; equations = {a*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(v[
        t, \ \[Theta], \ \[Phi]]\)\) + v[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(v[
        t, \ \[Theta], \ \[Phi]]\)\) + w[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(v[
        t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] == 
    w[t, \[Theta], \[Phi]]^2*Cos[\[Theta]]/Sin[\[Theta]] + 
     2*\[CapitalOmega]0*a0*w[t, \[Theta], \[Phi]]*
      Cos[\[Theta]] + (\[CapitalOmega]0*
        a0)^2 Sin[\[Theta]] Cos[\[Theta]] - g*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(h[
        t, \ \[Theta], \ \[Phi]]\)\) - \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(patmFullFunction[
        t, \ \[Theta], \ \[Phi]]\)\)*1/\[Rho], a*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(w[
        t, \ \[Theta], \ \[Phi]]\)\) + v[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(w[
        t, \ \[Theta], \ \[Phi]]\)\) + w[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(w[
        t, \ \[Theta], \ \[Phi]]\)\)/
      Sin[\[Theta]] == -w[t, \[Theta], \[Phi]]*v[t, \[Theta], \[Phi]]*
      Cos[\[Theta]]/Sin[\[Theta]] - 
     2*\[CapitalOmega]0*a0*v[t, \[Theta], \[Phi]]*Cos[\[Theta]] - 
     g*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(h[
        t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] - \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(patmFullFunction[
        t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]]*1/\[Rho], a*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(h[
        t, \ \[Theta], \ \[Phi]]\)\) + v[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(h[
        t, \ \[Theta], \ \[Phi]]\)\) + w[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(h[
        t, \ \[Theta], \ \[Phi]]\)\)/
      Sin[\[Theta]] + (h[t, \[Theta], \[Phi]] + H0 - 
        profileFullFunction[t, \[Theta], \[Phi]])*(\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(w[
          t, \ \[Theta], \ \[Phi]]\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\((v[
            t, \ \[Theta], \ \[Phi]]*Sin[\[Theta]])\)\))/
      Sin[\[Theta]] - v[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Theta]\)]\(profileFullFunction[
        t, \ \[Theta], \ \[Phi]]\)\) - w[t, \[Theta], \[Phi]]*\!\(
\*SubscriptBox[\(\[PartialD]\), \(\[Phi]\)]\(profileFullFunction[
        t, \ \[Theta], \ \[Phi]]\)\)/Sin[\[Theta]] == a*\!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]\(profileFullFunction[
       t, \ \[Theta], \ \[Phi]]\)\)} // Quiet;

initialconditions = {v[0, \[Theta], \[Phi]] == 
    v0Function[\[Theta], \[Phi]], 
   w[0, \[Theta], \[Phi]] == w0Function[\[Theta], \[Phi]], 
   h[0, \[Theta], \[Phi]] == h0Function[\[Theta], \[Phi]]};

eps = 10^-6; bc = {h[t, eps, \[Phi]] == 0, v[t, eps, \[Phi]] == 0, 
  w[t, eps, \[Phi]] == 0, h[t, Pi - eps, \[Phi]] == 0, 
  v[t, Pi - eps, \[Phi]] == 0, w[t, Pi - eps, \[Phi]] == 0, 
  v[t, \[Theta], 0] == v[t, \[Theta], 2 Pi], 
  w[t, \[Theta], 0] == w[t, \[Theta], 2 Pi], 
  h[t, \[Theta], 0] == h[t, \[Theta], 2 Pi]};

mol[n_Integer, o_ : "Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

sol = NDSolveValue[{equations, initialconditions, bc}, {v, w, h}, {t, 
    0, 1}, {\[Theta], eps, Pi - eps}, { \[Phi], 0, 2 Pi} , 
   Method -> mol[81, 2]] // AbsoluteTiming

It takes about 2180s on my Silver Pentium. Visualization

framex = 
  Table[Plot[
    sol[[2, 3]][t, x, Pi] - h0Function[x, y], {x, 10^-6, 3.14}, 
    PlotRange -> {-40, 40}, Axes -> False, Frame -> True], {t, 
    0.01, .3, 10^-2}];

framey = 
      Table[Plot[
        sol[[2, 3]][t, Pi/2, y] - h0Function[x, y], {y, 0, 2  Pi}, 
        PlotRange -> {-40, 40}, Axes -> False, Frame -> True], {t, 
        0.01, .3, 10^-2}];

In 2D it looks very funny

frames = 
 Table[DensityPlot[
   sol[[2, 3]][t, x, y] - h0Function[x, y], {y, 0, 2  Pi}, {x, 10^-6, 
    3.14}, PlotRange -> {-40, 40}, ColorFunction -> Hue, 
   Frame -> False, AspectRatio -> Automatic, MaxRecursion -> 2, 
   PlotPoints -> 20], {t, 0.03, .3, 10^-2}]

Answer 2

Its difficult or hopeless to NDSolve waves on a sphere with boundaries because the governing Laplacian has discrete spectrum. Eigenfunctions of the Laplacian with non-integer eigenvalues belong to parts of the sphere with haed boundaries along some coordinate lines.

The latitude periodic boundary for $\phi$ condition is quite harmless, the equation at poles requires a power factor of $\sin ^m \frac{\theta}{2}$ in order to generate a differentiable smooth function $\sin(m \phi) \ \sin^m\left(\frac{\theta}{2}\right)$.

In a wave equation, any discetization produces errors, that propagate and are tsunamis, too, as you observe in your simulation.

So it's of fundamental importance, to work inside the space of spherical harmonics, representing the discrete solution space in the same way as one uses Fourier methods on a torus: Fourier decompose, time evolute, Fourier compose.