sandbox/lbattershill/multilayer-slide/slide_layered.c

    Tsunami generation by a “pyroclastic” density current

    This is a simulation of a pyroclastic density current entering water and generating tsunami waves. We simulate a dense layer at base of flow, overlayed by a lighter layer. The setup is dimensionless and based on the experimental geometry of Bougouin et al., 2020, where a laboratory fluidised granular flow runs down a ramp and into water. This is the multilayer version of the Navier-Stokes case (Battershill et al., 2021).

    The setup is adapted from overflow.c to simulate a flow of variable density layers into water, using the Boussinesq buoyancy module of the non-hydrostatic layered scheme.

    Closeup view: the colorscale is the relative density variation.

    Full tank view: the colorscale is the relative density variation.

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

    Setting \Delta \rho directly.

    #define drho(T) T

    Including the relevant modules and intialising the simulation

    #include "layered/dr.h"
    
    #include "layered/remap.h"
    #include "layered/rpe.h"
    #include "layered/perfs.h"
    
    #define nlayers 26
    
    double rho_0 = 0.1;
    double rho_2 = -0.1;
    double rho_3 = 0.;
    
    //double nu_H = 0.0001; // 1e-2;
    
    int main()
    {
      size (32);
      //breaking = 0.7;
      N = 512; //2048
      nl = nlayers;
      nu = 0.0001;
      DT = 0.0025;
      G = 1.;
      cell_lim = mono_limit;
      system ("rm -r gnuplot");
      system ("mkdir gnuplot");
      system ("mkdir gnuplot/all");
      system ("mkdir gnuplot/closeup");
      system ("mkdir gnuplot/closeup/density");
      system ("mkdir gnuplot/closeup/ux");
      const scalar slip[] = HUGE;
    incompatible types when assigning to type ‘vector’ {aka ‘struct ’} from type ‘scalar’ {aka ‘const struct ’}
      lambda_b = slip;
      run();
    }

    Define the geometry.

    #define theta0 0.268 //angle
    #define width 0.338/0.265 //width of flow
    #define gradient_m -1.*tan(theta0) //gradient of slope
    #define Hi 1. //initial water depth
    #define slope_exposed 1.
    #define Hel (slope_exposed/0.265)*sin(theta0)
    #define c_intercept tan(theta0)*(width+ (Hi+Hel)/tan(theta0)) - Hi
    #define height_cm 0.185/0.265
    #define frac 6./7.
    #define c_intercept2 height_cm + ((frac*height_cm)/(1- frac))

    Initialise domain.

    event init (i = 0)
    {
      foreach() {
          zb[] =  x < width ? Hel : (x  < width + (Hi+Hel)/tan(theta0) ?
          			   gradient_m*x + c_intercept  : -Hi); 
        foreach_layer() {
          if (x <= width*frac) {// Granular fluid
    	// DEPTH
    	h[] =  (height_cm)/(nl);
    	// BUOYANCY
    	if (point.l < nl/2) 
    	  T[] = rho_0;
    	else
    	  T[] = rho_2;
    	
          }
          else if (x < width) { //Granular fluid
    	h[] = (-1*x*(height_cm/(width-(width*frac))) + c_intercept2 ) /nl;	
    	if (point.l < nl/2)
    	  T[] = rho_0;
    	else
    	  T[] = rho_2;
    	
          }
          else if (x < width + Hel/tan(theta0)) { // Gap
    	// DEPTH
    	h[] = 0.00002;
    	// BUOYANCY
    	T[] = rho_0;
          }
          else { //Water
    	// DEPTH
    	h[] = - zb[]/nl;
    	// BUOYANCY
    	T[] = rho_3;
          }
        }
      }
    }

    Outputs

    event logfile (i += 10)
    {
      static double rpe0 = 0., rpen = 0., tn = 0.;
      double rpe = RPE();
      double PE, KE;
      energy (&PE, &KE);
      if (i == 0) {
        rpe0 = rpe;
        fprintf (stderr, "t  dt   rpe-rpe0   rpe0   PE   KE   d_t(rpe)   mgp.i\n");
      }
      fprintf (stderr, "%g %g %.12g %.12g %.12g %.12g %.12g %d\n", t, dt,
    	   rpe - rpe0, rpe0, PE, KE,
    	   t > tn ? (rpe - rpen)/(t - tn)/L0 : 0.,
    
    #if NH	   
    	   mgp.i
    #else
    	   mgH.i
    #endif
    	   );
      rpen = rpe, tn = t;
    }

    Animations including density, u.x…

    void setup (FILE * fp)
    {
      fprintf (fp,
    #if ISOPYCNAL
    	   "set pm3d map corners2color c2\n"
    #else
    	   "set pm3d map\n"
    #endif
    	   "# jet colormap\n"
    	   "set palette defined ( 0 0 0 0.5647, 0.125 0 0.05882 1, 0.25"
    	   " 0 0.5647 1, 0.375 0.05882 1 0.9333, 0.5 0.5647 1 0.4392, 0.625"
    	   " 1 0.9333 0, 0.75 1 0.4392 0, 0.875 0.9333 0 0, 1 0.498 0 0 )\n"
    	   "unset key\n"
    	   "set xlabel 'x'\n"
    	   "set ylabel 'depth'\n"
    	   );
    }
    
    void plot_ux (FILE * fp)
    {
      fprintf (fp,
    	   "set cbrange [-0.1:30]\n" // u.x
    	   "set title 't = %.2f'\n"
    	   "sp '-' u 1:2:3\n",t);
      foreach (serial) {
        double z = zb[];
        fprintf (fp, "%g %g %g %g\n", x, z, u.x[], T[]);
        foreach_layer() {
          z += h[];
          fprintf (fp, "%g %g %g %g\n", x, z, u.x[], T[]);
        }
        fprintf (fp, "\n");
      }
      fprintf (fp, "e\n\n");
      fflush (fp);
    }
    
    void plot_density (FILE * fp)
    {
      fprintf (fp,
    	   "set cbrange [-0.1:0.1]\n" // Density
    	   "set title 't/T0 = %.2f'\n"
    	   "sp '-' u 1:2:4\n",t);
      foreach (serial) {
        double z = zb[];
        fprintf (fp, "%g %g %g %g\n", x, z, u.x[], T[]);
        foreach_layer() {
          z += h[];
          fprintf (fp, "%g %g %g %g\n", x, z, u.x[], T[]);
        }
        fprintf (fp, "\n");
      }
      fprintf (fp, "e\n\n");
      fflush (fp);
    }
    
    event gnuplot (t += 0.1)
    {
      static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
      if (i == 0)
        setup (fp);
      // Uncomment for animation during run.
      /*
        if (getenv ("DISPLAY")) {
        fprintf (fp, "set term x11\n");
        plot_density (fp);
        }
      */
      fprintf (fp,
    	   "set term pngcairo font \",10\" size 1600,620\n"
    	   "set xrange [0:12.5]\n"
    	   "set yrange [-1.5:4]\n"
    	   "set size ratio -1\n"
    	   "set ytics\n"
    	   "set output 'gnuplot/closeup/density/plot-%06d.png'\n", i);
      plot_density (fp);
      
      fprintf (fp,
    	   "set term pngcairo font \",10\" size 1600,266\n"
    	   "set xrange [0:32]\n"
    	   "set yrange [-1.5:2.5]\n"
    	   "set size ratio -1\n"
    	   "unset ytics\n"
    	   "set output 'gnuplot/all/plot-%06d.png'\n", i);
      plot_density (fp);  
    }

    Equally placed wave gauges.

    Gauge gauges[] = {
      {"WG_3m",  3./0.265, 0.},
      {"WG_4m", 4./0.265, 0.},
      {"WG_5m", 5./0.265, 0.},
      {"WG_6m", 6./0.265, 0.},
      {"WG_7m", 7./0.265, 0.},
      {NULL}
    };
    
    event output (t += 0.1)
      output_gauges (gauges, {eta});

    We plot the wave gauges:

    set xlabel 't/T0'
    set ylabel 'y'
    plot 'WG_3m' u 1:2 w l t 'x = 11.3','WG_4m' u 1:2 w l t 'x = 15.1', 'WG_5m' u 1:2 w l t 'x = 18.9'
    Evolution of the free surface elevation at x = (script)

    Evolution of the free surface elevation at x = (script)

    event end (t = 25.)
    { 
      system ("for f in gnuplot/all/plot-*.png; do"
    	  " convert $f ppm:- && rm -f $f; done | "
    	  "ppm2mp4 movie_normal.mp4");
      system ("for f in gnuplot/closeup/density/plot-*.png; do"
    	  " convert $f ppm:- && rm -f $f; done | "
    	  "ppm2mp4 movie_closeup.mp4");
      fprintf (stderr, "\n\nDone\n");
    }

    References

    [battershill2021]

    Lily Battershill, Colin Whittaker, Emily Lane, Stephane Popinet, James White, William Power, and P Nomikou. Numerical simulations of a fluidized granular flow entry into water: insights into modeling tsunami generation by pyroclastic density currents. 2021.

    [bougouin2020]

    Alexis Bougouin, Raphael Paris, and Olivier Roche. Impact of fluidized granular flows into water: implications for tsunamis generated by pyroclastic flows. Journal of Geophysical Research: Solid Earth, 125(5):e2019JB018954, 2020. [ DOI ]