src/test/large.c

    Large-amplitude standing wave

    This test case is discussed in Popinet (2020) and compares the results obtained with the layered model with the Navier-Stokes/VOF solver for the large-amplitude, non-linear oscillation of a gravity wave.

    set xlabel 'x'
    set ylabel 'z'
    plot 'log' u 1:2 w l t 'multilayer', \
         '../large-ns/log' u 1:2 w l t 'VOF', \
         'large.nometric' u 1:2 w l t 'multilayer (no metric)'
    Free surface profiles (script)

    Free surface profiles (script)

    set xlabel 't'
    set ylabel ''
    set yr [:0.8]
    plot 'out' u 1:2 w l t 'Maximum vertical velocity', \
         'out' u 1:3 w l t 'Maximum slope', \
         '../large-ns/out' u 1:2 every 10 t '', \
         '../large-ns/out' u 1:3 every 10 t '', \
         9.81*x t '9.81 t'
    Maximum vertical velocity and slope (script)

    Maximum vertical velocity and slope (script)

    See also

    #include "grid/multigrid1D.h"
    #include "layered/hydro.h"
    #include "layered/nh.h"
    #include "layered/remap.h"
    
    int main()
    {
      N = 512;
      G = 9.81;
      nl = 20;
      max_slope = 1.; // a bit better than with the more dissipative default 0.577
      CFL_H = 0.5;
      NITERMIN = 2;
      run();
    }
    
    event init (i = 0)
    {
      double k = 2.*pi, a = 0.07;
      foreach() {
        zb[] = - 0.5;
        foreach_layer()
          h[] = (a*cos(k*x) - zb[])/nl;
      }
    }
    
    #if 0
    event gnuplot (i += 1) {
      static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
      if (i == 0)
        fprintf (fp, "set term x11; set size ratio -1\n");
      fprintf (fp,
    	   "set title 'nl = %d, t = %.2f'\n"
    	   "p [%g:%g][-0.5:]'-' u 1:3:2 w filledcu lc 3 t '',"
    	   " '' u 1:(-1):3 t '' w filledcu lc -1", nl, t,
    	   X0, X0 + L0);
      int i = 4;
      foreach_layer()
        fprintf (fp, ", '' u 1:%d w l lw 2 t ''", i++);
      fprintf (fp, "\n");
      foreach (serial) {
        double H = 0.;
        foreach_layer()
          H += h[];
        fprintf (fp, "%g %g %g", x, zb[] + H, zb[]);
        double z = zb[];
        foreach_layer() {
          fprintf (fp, " %g", z);
          z += h[];
        }
        fprintf (fp, "\n");
      }
      fprintf (fp, "e\n\n");
      //  fprintf (fp, "pause 0.01\n");
      fflush (fp);
    }
    #endif
    
    event logfile (i++)
    {
      double wmax = 0., smax = 0.[0];
      foreach (reduction (max:wmax) reduction (max:smax)) {
        double Hm = 0., Hp = 0.;
        foreach_layer() {
          if (w[] > wmax)
    	wmax = w[];
          Hm += h[-1], Hp += h[1];
        }
        if ((Hp - Hm)/(2.*Delta) > smax)
          smax = (Hp - Hm)/(2.*Delta);
      }
      printf ("%g %g %g\n", t, wmax, smax);
    }
    
    event profiles (t = 0.1; t += 0.1; t <= 0.5)
    {
      foreach()
        fprintf (stderr, "%g %g\n", x, eta[]);
      fprintf (stderr, "\n");
    }