src/test/galilean_invariance.c

    Advection of a rippled interface

    We test the ability of the multilayer solver to transport a corrugated interface with a constant velocity profile, i.e. we check that the solver preserves the property of Galilean invariance.

    Advection of a rippled interface

    Advection of a rippled interface

    In the absence of body forces, the interface is initialised with a sinusoidal shape which is then swept away with a uniform velocity profile U (this moving interface is therefore not to be confused with a inertial-gravity propagating wave).

    The resulting solid-body like motion is an exact solution of the multilayer set of equations with the non-hydrostatic corrections, with \phi = 0.

    #include "grid/multigrid1D.h"
    #include "layered/hydro.h"
    #include "layered/nh.h"
    #include "layered/remap.h"

    We non-dimensionalise the problem with the advection velocity U and the wavelength \lambda, and we set the amplitude of the wave to one tenth of the channel depth.

    #define wavelength 1.
    #define U 1.
    #define wavenumber (2.*pi/wavelength)
    #define depth 0.5
    #define amplitude (depth/10.)
    #define endoftime (1.*wavelength/U)
    
    double maxwaveerror;

    Main function

    The boundary conditions are set to periodic. The test is done in a zero-gravity environment. The default minmod slope limiter is turned off to avoid blunting off the wave crests.

    int main()
    {
      periodic (right);
      G = 0.;
      gradient = NULL;
    
      for (N = 64; N <= 256; N *= 2)
        for (nl = 1; nl <= 16; nl *= 2)
          run();
    }

    Initialisation

    The initial shape of the interface is a simple sine, and the velocity field is set to a constant velocity U = 1 and w = 0 everywhere.

    event init (i = 0)
    {
      foreach() {
        double H = depth + amplitude*cos(wavenumber*x);
        foreach_layer() {
          h[] = H/nl;
          u.x[] = U;
        }
      }
      maxwaveerror = 0.;
    }

    Outputs

    We keep track of the amplitude of the wave vs time, and for a typical case (N = 128 and \text{nl} = 4) we export deeper checks on the non-hydrostatic pressure \phi and velocity u.

    event monitoring (t += endoftime/100; t = 0.; t <= endoftime)
    {
      stats s = statsf(eta);
      double etaamp = fabs((s.max - s.min)/(2*amplitude) - 1);
      if (N == 128 && nl == 4) {
        char name[80];
        sprintf (name, "output-N-%d-nl-%d", N, nl);
        static FILE * fpout = fopen (name, "w");
    
        scalar absdevU[];
        scalar absPhi[];
        foreach()
          foreach_layer() {
            absdevU[] = fabs(u.x[]-U);
            absPhi[]  = fabs(phi[]);
          }
        boundary({absdevU,absPhi});
        double maxdevU = statsf(absdevU).max/U;
        double maxPhi = statsf(absPhi).max;
        assert (maxPhi < 1e-14 && maxdevU < 1e-14);
        fprintf (fpout, "%g %g %g %g\n", t, maxPhi, maxdevU, etaamp);
      }
      maxwaveerror = max(maxwaveerror, etaamp);
    }

    The reference file of this test is just based on the overall maximum wave amplitude error.

    event errorlog (t = end)
    {
      fprintf (stderr, "%d %d %.3e\n", N, nl, maxwaveerror);
    }

    Movie

    This is how we export the headline animated gif.

    event movie (i += 5; t <= (wavelength/U))
    {
      if (N == 128 && nl == 4) {
        static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
        if (i == 0)
          fprintf (fp,
    	       "set term pngcairo enhanced font ',10' size 800,250\n"
    	       "set size ratio -1\n"
    	       "unset key\n");
        fprintf (fp,
    	     "set output 'plot%04d.png'\n"
    	     "set title 't = %.2f'\n"	     
    	     "p [%g:%g][0:]'-' u 1:3:2 w filledcu lt rgb \"#993498DB\", "
    	     " '' u 1:2 w l lw 4 lt rgb \"#3498DB\" ",
    	     i/5, t, X0, X0 + L0);
        int i = 4;
        foreach_layer()
          fprintf (fp, ", '' u 1:%d w l dt 2 lt rgb \"gray50\"", i++);
        fprintf (fp, "\n");
        foreach_leaf() {
          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");
        fflush (fp);
      }
    }
    
    event moviemaker (t = end)
    {
      if (N == 128 && nl == 4)
        system ("mogrify -format gif plot*.png && "
          "gifsicle --delay 10 --colors 256 --loop plot*.gif > wave.gif && "
          "rm -f plot*.*");
    }

    Results

    For the typical case N = 128 and \text{nl} = 4 the time history of the Poisson error \|\phi\|_\infty and of any deviation of the velocity to the setpoint \|u-U\|_\infty – which would be responsible for phase errors – are represented and are observed to be of the order of the machine precision.

    set term @SVG size 960,240 font ',10'
    set size ratio 2./(1. + sqrt(5.))
    set logscale y
    set grid
    set style line 1 \
        linecolor rgb '#774F38' \
        linetype 1 linewidth 2 \
        pointtype 7 pointsize 1.
    set style line 2 \
        linecolor rgb '#E08E79' \
        linetype 1 linewidth 2 \
        pointtype 7 pointsize 1.
    set style line 3 \
        linecolor rgb '#F1D4AF' \
        linetype 1 linewidth 2 \
        pointtype 7 pointsize 1.
    set style line 4 \
        linecolor rgb '#ECE5CE' \
        linetype 1 linewidth 2 \
        pointtype 7 pointsize 1.
    set style line 5 \
        linecolor rgb '#C5E0DC' \
        linetype 1 linewidth 2 \
        pointtype 7 pointsize 1.
    
    set ylabel 'Error (phi)'
    set xlabel "t"
    set key left
    set multiplot layout 1, 2
    plot [0:1][1e-16:1e-10]'./output-N-128-nl-4' u 1:2 t 'Error on phi' with lp ls 1
    set ylabel 'Error (u)'
    plot [0:1][1e-16:1e-10]'./output-N-128-nl-4' u 1:3 t 'Error on u' with lp ls 2
    unset multiplot
    Left: error on \phi vs time estimated as \|\phi\|_\infty (Poisson error). Right: error on horizontal velocity \|u-U\|_\infty vs time (phase error) (script)

    Left: error on \phi vs time estimated as \|\phi\|_\infty (Poisson error). Right: error on horizontal velocity \|u-U\|_\infty vs time (phase error) (script)

    For the same case (N = 128 and \text{nl} = 4) the relative error in wave amplitude (dissipation error) is monitored through time, and is seen to be of the order of 10^{-4}.

    set term @SVG size 640,320 font ',10'
    set ylabel 'Error (amplitude)'
    plot './output-N-128-nl-4' u 1:4 t 'Relative error on amplitude' with lp ls 3
    Relative error of the wave amplitude with time (dissipation error). (script)

    Relative error of the wave amplitude with time (dissipation error). (script)

    Convergence

    We track the value of the relative error on wave amplitude for various number of gridpoints N and number of layers \text{nl}. For up to 16 stacked layers the results are undistinguishable, and the error is observed to decrease as N^{-2}, indicating a second-order precision.

    set xlabel 'Number of grid points'
    set ylabel 'Error'
    set key bottom left
    set logscale x 2
    set label 1 "N^{-2}" at 180,5e1*180**-2 font ',12' textcolor rgb '#E08E79' offset 0,1
    filter1 = '< awk ''{ if ($2 == 1) print $1, $2, $3}'' '
    filter2 = '< awk ''{ if ($2 == 2) print $1, $2, $3}'' '
    filter4 = '< awk ''{ if ($2 == 4) print $1, $2, $3}'' '
    filter8 = '< awk ''{ if ($2 == 8) print $1, $2, $3}'' '
    filter16 = '< awk ''{ if ($2 == 16) print $1, $2, $3}'' '
    plot filter1.'log' u 1:3 t "nl = 1" with linespoints linestyle 1, \
         filter2.'log' u 1:3 t "nl = 2" with linespoints linestyle 2, \
         filter4.'log' u 1:3 t "nl = 4" with linespoints linestyle 3, \
         filter8.'log' u 1:3 t "nl = 8" with linespoints linestyle 4, \
         filter16.'log' u 1:3 t "nl = 16" with linespoints linestyle 5, \
         [96:96<<2]5e1*x**-2 t '' with lines linestyle 2 lw 3
    Variation of the relative error on wave amplitude with horizontal resolution, for different (fixed) number of layers. (script)

    Variation of the relative error on wave amplitude with horizontal resolution, for different (fixed) number of layers. (script)