/** # Solitary wave run-up on a plane beach We use the [Green-Naghdi](/src/green-naghdi.h) or the [layered](/src/layered/hydro.h) solver to reproduce this classical test case based on the experiments of [Synolakis, 1987](/src/references.bib#synolakis1987). */ #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](/src/references.bib#yamazaki2009)). 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](soliton.c) 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](/src/references.bib#yamazaki2009) 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](/src/references.bib#yamazaki2009) Figure 6, obtained with a different set of depth-averaged equations and the results of [Bonneton et al, 2011](/src/references.bib#bonneton2011), 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). ~~~gnuplot Comparison of model predictions and experimental snapshots 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 ~~~ */