# src/examples/shoal.c

# Periodic wave propagation over an ellipsoidal shoal

We follow Lannes and Marche, 2014 and try to reproduce the experiment of Berkhoff et al, 1982. The numerical wave tank is 25^{2} metres and periodic waves are generated on the left-hand side and are damped on the right-hand side.

```
#include "grid/multigrid.h"
#include "green-naghdi.h"
int main()
{
X0 = -10;
L0 = 25;
Y0 = -L0/2.;
G = 9.81;
N = 1024;
```

We turn off limiting to try to reduce wave dissipation.

```
gradient = NULL;
run();
}
```

Periodic waves with period one second are generated on the left-hand-side. We tune the amplitude of the “radiation” condition to match that of the experiment as measured by wave gauges.

```
u.n[left] = - radiation (0.042*sin(2.*pi*t/1.));
u.n[right] = + radiation (0);
```

We declare a new field to store the maximum wave amplitude.

```
scalar maxa[];
event init (i = 0)
{
```

The bathymetry is an inclined and skewed plane combined with an ellipsoidal shoal.

```
double h0 = 0.45;
double cosa = cos (20.*π/180.), sina = sin (20.*π/180.);
foreach() {
double xr = x*cosa - y*sina, yr = x*sina + y*cosa;
double z0 = xr >= -5.82 ? (5.82 + xr)/50. : 0.;
double zs = sq(xr/3.) + sq(yr/4.) <= 1. ?
-0.3 + 0.5*sqrt(1. - sq(xr/3.75) - sq(yr/5.)) : 0.;
zb[] = z0 + zs - h0;
h[] = max(-zb[], 0.);
maxa[] = 0.;
}
}
```

To implement an absorbing boundary condition, we add an area for $x>12$ for which quadratic friction increases linearly with $x$.

```
event friction (i++) {
foreach()
if (x > 12.) {
double a = h[] < dry ? HUGE : 1. + 2.*(x - 12.)*dt*norm(u)/h[];
foreach_dimension()
u.x[] /= a;
}
}
```

Optionally, we can make a “stroboscopic” movie of the wave field. This is useful to check the amount of waves reflected from the outflow.

```
#if 0
event movie (t += 1) {
output_ppm (η, min = -0.04, max = 0.04, n = 512);
}
#endif
```

After the wave field is established ($t>=40$) we record the maximum wave amplitude.

```
event maximum (t = 40; i++) {
foreach()
if (fabs(η[]) > maxa[])
maxa[] = fabs(η[]);
}
```

At the end of the simulation, we output the maximum wave amplitudes along the cross-sections corresponding with the experimental data.

```
event end (t = 50) {
FILE * fp = fopen ("end", "w");
output_field ({η,zb,maxa}, fp, linear = true);
FILE * fp2 = fopen ("section2", "w");
FILE * fp5 = fopen ("section5", "w");
FILE * fp7 = fopen ("section7", "w");
for (double y = -10.; y <= 10.; y += 0.02) {
fprintf (fp2, "%g %g\n", y, interpolate (maxa, 3., y));
fprintf (stderr, "%g %g\n", y, interpolate (maxa, 5., y));
fprintf (fp5, "%g %g\n", y, interpolate (maxa, 9., y));
}
for (double x = -10.; x <= 10.; x += 0.02)
fprintf (fp7, "%g %g\n", x, interpolate (maxa, x, 0.));
}
```

The results are comparable to Lannes and Marche, 2014, Figure 19, but we have to use a higher resolution (1024^{2} instead of 300^{2}) because our numerical scheme is only second order (versus 4th order).