src/test/tsunami.c

    Typical (1D) tsunami wave

    The initial surface displacement is a Gaussian bump with an amplitude of 10 metres and caracteristic width 25 km. The water depth is constant at 7000 metres. These parameters are fairly typical of fault-triggered tsunamis close to subduction zones.

    As illustrated below, dispersive effects are important.

    Note also that older versions of the non-hydrostatic multilayer solver were affected by a grid-scale instability for this test case.

    set xlabel 'x (km)'
    set ylabel 'Amplitude (m)'
    unset key
    plot 'log' u ($1/1e3):2 w l
    Wave shape after one hour (script)

    Wave shape after one hour (script)

    set ylabel 'Horizontal velocity (m/s)'
    plot 'log' u ($1/1e3):3 w l
    Horizontal velocity after one hour (script)

    Horizontal velocity after one hour (script)

    set ylabel 'Vertical velocity (m/s)'
    plot 'log' u ($1/1e3):4 w l
    Vertical velocity after one hour (script)

    Vertical velocity after one hour (script)

    #include "grid/multigrid1D.h"
    #include "layered/hydro.h"
    #include "layered/nh.h"
    #include "layered/check_eta.h"
    
    int main()
    {
      size (2500e3);
      origin (-L0/2.);
      G = 9.81;
      N = 2 << 12;
      CFL_H = 0.5;
      run();
    }
    
    event init (i = 0)
    {
      res_eta = new scalar;
      conserve_elevation();
      foreach() {
        zb[] = - 7000;
        h[] = - zb[] + 10.*exp(-sq(x/25e3));
      }
      boundary ({h});
    }
    
    #if 1
    event gnuplot (i += 10)
    {
      static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
      if (i == 0)
        fprintf (fp, "set term x11\n"
    	     "set xlabel 'x (km)'\n"
    	     "set ylabel 'Wave height (m)'\n"
    	     //	     "set yrange [-4:10]\n"
    	     "set xrange [%g:%g]\n"
    	     // "unset key\n"
    	     , X0/1e3, (X0 + L0)/1e3);
      FILE * fp1 = fopen ("gnuplot", "w");
      foreach()
        fprintf (fp1, "%g %g %g %g\n", x, eta[], deta[], res_eta[]);
      fclose (fp1);
      fprintf (fp,
    	   "set title 't = %.2f min'\n"
    	   "p 'gnuplot' u ($1/1e3):3 w l t 'etap - eta',"
    	   "  '' u ($1/1e3):4 w d t 'res'\n"
    	   //	   "pause 1\n"
    	   , t/60.);
      fflush (fp);
    }
    #endif
    
    event plots (t = 3600)
    {
      foreach()
        fprintf (stderr, "%g %.4g %.4g %.4g\n", x, eta[], u.x[], w[]);
      fprintf (stderr, "\n");
    }