src/test/beach-ml.c

    Solitary wave run-up on a plane beach

    We use the Green-Naghdi or the layered solver to reproduce this classical test case based on the experiments of Synolakis, 1987.

    #include "grid/multigrid1D.h"
#if ML
#  include "layered/hydro.h"
#  include "layered/nh.h"
#  include "layered/remap.h"
#else
#  include "green-naghdi.h"
#endif

    The problem is non-dimensionalised by the water depth h_0 and the acceleration of gravity. The amplitude of the solitary wave is 0.28 and the slope of the beach is 1/19.85. The length L is the estimated wavelength of the solitary wave (see section 4.3 of Yamazaki et al, 2009). We set the coordinate system (X0 and L0) so that the origin is the intersection of the beach with the water level.

    double h0 = 1., A = 0.28, L;
double slope = 1./19.85;

int main()
{
  double k = sqrt(3.*A/4/cube(h0));
  L = 2./k*acosh(sqrt(1./0.05));
  X0 = - h0/slope - L/2. - L;
  L0 = 6.*L;
  N = 1024;
  G = 1.;
#if ML
  nl = 2;
  breaking = 0.07;
#else
  alpha_d = 1.;

    We try to tune the “breaking slope” to get better agreement for t=20.

      breaking = 0.4;
#endif
  run();  
}

    The initial wave is the analytical soliton for the Green-Naghdi equations.

    double sech2 (double x) {
  double a = 2./(exp(x) + exp(-x));
  return a*a;
}

double soliton (double x, double t)
{
  double c = sqrt(G*(1. + A)*h0), psi = x - c*t;
  double k = sqrt(3.*A*h0)/(2.*h0*sqrt(h0*(1. + A)));
  return A*h0*sech2 (k*psi);
}

    See figure 5 of Yamazaki et al, 2009 for the definition of the bathymetry.

    event init (i = 0)
{
  double c = sqrt(G*(1. + A)*h0);
  foreach() {
    double eta = soliton (x + h0/slope + L/2., t);
    zb[] = max (slope*x, -h0);
#if ML
    foreach_layer() {
      h[] = max (0., eta - zb[])/nl;
      u.x[] = c*eta/(h0 + eta);
    }
#else
    h[] = max (0., eta - zb[]);
    u.x[] = c*eta/(h0 + eta);
#endif
  }
}

    Friction is important for this test case. We implement a simple time-implicit quadratic bottom friction and tune the coefficient to obtain a runup comparable with the experiment.

    event friction (i++) {
  foreach() {
#if ML
    double Q = 0., H = 0.;
    foreach_layer()
      H += h[], Q += h[]*u.x[];
    double a = H < dry ? HUGE : 1. + 5e-3*dt*fabs(Q)/sq(H);
    foreach_layer()
      u.x[] /= a;
#else
    double a = h[] < dry ? HUGE : 1. + 5e-3*dt*norm(u)/h[];
    foreach_dimension()
      u.x[] /= a;
#endif
  }
}

    We use gnuplot to display an animation while the simulation is running.

    event gnuplot (i += 5) {
  static FILE * fp = popen ("gnuplot 2> /dev/null", "w");
  fprintf (fp,
	   "set title 't = %.2f'\n"
	   "p [-20:12][-0.2:0.6]'-' u 1:3:2 w filledcu lc 3 t '',"
	   " '' u 1:(-1):3 t '' w filledcu lc -1\n", t);
  foreach()    
    fprintf (fp, "%g %g %g\n", x, eta[], zb[]);
  fprintf (fp, "e\n\n");
  fprintf (stderr, "%.3f %.3f\n", t, statsf(u.x).max);
}

    We output profiles at the same times as the experimental data, in separate files indexed by the time.

    event output (t <= 65; t += 5) {
  char name[80];
  sprintf (name, "out-%g", t);
  FILE * fp = fopen (name, "w");
  foreach()
    fprintf (fp, "%g %g %g\n", x, eta[], zb[]);
  fprintf (fp, "\n");
  fclose (fp);
}

    The agreement with the experiments (circles) is satisfactory and is comparable to the numerical results of Yamazaki et al, 2009 Figure 6, obtained with a different set of depth-averaged equations and the results of Bonneton et al, 2011, Figure 8, although this latter model seems to do a better job of capturing breaking around t=20 (they use a more sophisticated breaking criterion).

    set term svg enhanced size 640,640 font ",10"
set xrange [-20:12]
set yrange [-0.2:0.6]
set ytics -0.2,0.2,0.6
set key top left
set multiplot layout 6,2 scale 1.05,1.1
set rmargin 2
set tmargin 0.5
unset xtics
plot 'out-10' w l lw 2 t 't=10', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-10' pt 6 lc -1 ps 0.5 t ''
unset ytics
plot 'out-15' w l lw 2 t 't=15', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-15' pt 6 lc -1 ps 0.5 t ''
set ytics -0.2,0.2,0.6
plot 'out-20' w l lw 2 t 't=20', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-20' pt 6 lc -1 ps 0.5 t ''
unset ytics
plot 'out-25' w l lw 2 t 't=25', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-25' pt 6 lc -1 ps 0.5 t ''
set ytics -0.2,0.2,0.6
plot 'out-30' w l lw 2 t 't=30', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-30' pt 6 lc -1 ps 0.5 t ''
unset ytics
plot 'out-35' w l lw 2 t 't=35', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-35' pt 6 lc -1 ps 0.5 t ''
set ytics -0.2,0.2,0.6
plot 'out-40' w l lw 2 t 't=40', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-40' pt 6 lc -1 ps 0.5 t ''
unset ytics
plot 'out-45' w l lw 2 t 't=45', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-45' pt 6 lc -1 ps 0.5 t ''
set ytics -0.2,0.2,0.6
plot 'out-50' w l lw 2 t 't=50', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-50' pt 6 lc -1 ps 0.5 t ''
unset ytics
plot 'out-55' w l lw 2 t 't=55', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-55' pt 6 lc -1 ps 0.5 t ''
set ytics -0.2,0.2,0.6
set xtics
plot 'out-60' w l lw 2 t 't=60', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-60' pt 6 lc -1 ps 0.5 t ''
unset ytics
plot 'out-65' w l lw 2 t 't=65', '' u 1:3 w filledcu x1 lc -1 t '', \
     '../t-65' pt 6 lc -1 ps 0.5 t ''
unset multiplot
    Comparison of model predictions and experimental snapshots (script)

    Comparison of model predictions and experimental snapshots (script)