# src/test/conical.c

# Runup of a solitary wave on a conical island

In this test case, we compare experimental data with solutions of the Saint-Venant and Green-Naghdi equations. For details on the test case and experimental data many references are available (e.g. Hou el al, 2013, Nikolos and Delis, 2009, Lannes and Marche, 2014, Kazolea et al, 2012 etc…).

We also test both the explicit and the implicit versions of the Saint-Venant solver.

```
#if SAINT_VENANT
# if IMPLICIT
# include "saint-venant-implicit.h"
# else
# include "saint-venant.h"
# endif
# define MGD 0
#else
# include "green-naghdi.h"
# define MGD mgD.i
#endif
```

The domain is 27.6 metres squared, resolved with a maximum of 256^{2} grid points.

```
#define MAXLEVEL 8
#define MINLEVEL 0
int main()
{
G = 9.81;
size (27.6);
```

We set an “acoustic” CFL for the implicit solver, so that the timestep is comparable to that of the explicit Saint-Venant solver. With larger timesteps the scheme is too dissipative.

```
#if IMPLICIT
CFLa = 0.25;
#endif
init_grid (1 << MAXLEVEL);
run();
}
```

We reproduce case B i.e. a relative amplitude of the solitary wave of 0.09.

```
double H0 = 0.32, H = 0.32*0.09, T = 2.45;
#define C sqrt(G*(H0 + H))
double sech2 (double x)
{
double e = exp (-x);
double s = 2.*e/(1. + e*e);
return s*s;
}
double eta0 (double t)
{
return H*sech2 (sqrt(3.*H/(4.*cube(H0)))*C*(t - T));
}
double uleft (double t)
{
double e = eta0 (t);
return C*e/(H0 + e);
}
```

We use the definition of the solitary wave to impose conditions on the left side (i.e. the “wave paddle”).

```
h[left] = H0 + eta0(t);
#if IMPLICIT
q.n[left] = (H0 + eta0(t))*uleft (t);
q.t[left] = 0.;
#else
u.n[left] = uleft (t);
u.t[left] = 0.;
#endif
```

This is the definition of the topography of the conical island.

```
double island (double x, double y)
{
double r0 = 3.6, r1 = 1.1, H0 = 0.625;
x -= 25.92/2.; y -= L0/2.;
double r = sqrt(x*x + y*y);
return r > r0 ? 0. : r < r1 ? H0 : (r - r0)/(r1 - r0)*H0;
}
```

On the tree we need to make sure that refined cells are initialised with the correct topography.

```
#if TREE
void refine_zb (Point point, scalar zb)
{
foreach_child()
zb[] = island (x, y);
}
#endif
event init (i = 0)
{
#if TREE
```

We setup the refinement function and choose to conserve surface elevation rather than volume.

```
zb.refine = refine_zb;
conserve_elevation();
#endif
```

These are the initial conditions i.e. conical island and lake-at-rest.

```
foreach() {
zb[] = island (x, y);
h[] = max(0., H0 - zb[]);
}
}
```

We define the wave gauges corresponding with the experimental locations.

```
Gauge gauges[] = {
{"WG3", 6.82, 13.05},
{"WG6", 9.36, 13.80},
{"WG9", 10.36, 13.80},
{"WG16", 12.96, 11.22},
{"WG22", 15.56, 13.80},
{NULL}
};
```

We output snapshots of fields at various times.

```
event outputfile (t = {9, 12, 13, 14, 20})
{
static int nf = 0;
printf ("file: conical-%d\n", nf);
output_field ({h,zb}, stdout, N, linear = true);
```

Here we output the level of refinement of the grid.

```
printf ("file: level-%d\n", nf);
scalar l[];
foreach()
l[] = level;
output_field ({l}, stdout, N);
nf++;
}
```

Which gives the following wave (left column) and mesh (right column) evolutions.

`event logfile (i++) {`

We output various diagnostics.

```
stats s = statsf (h);
#if IMPLICIT
scalar u[];
foreach()
u[] = h[] > dry ? q.x[]/h[] : 0.;
norm n = normf (u);
#else
norm n = normf (u.x);
#endif
if (i == 0)
fprintf (stderr, "t i h.min h.max h.sum u.x.rms u.x.max dt mgD.i\n");
fprintf (stderr, "%g %d %.8f %g %.8f %g %.4g %g %d\n",
t, i, s.min, s.max, s.sum, n.rms, n.max, dt, MGD);
```

Here we output surface elevation at the various gauges…

```
#if IMPLICIT
scalar η[];
foreach()
η[] = h[] > dry ? h[] + zb[] : nodata;
boundary ({η});
#endif
output_gauges (gauges, {η});
}
```

… and compare it with experimental data for the Green-Naghdi solver (blue) and the Saint-Venant solver (magenta). The results for the Green-Naghdi solver are very similar to those of Kazolea et al, 2012, Figure 14 and Lannes and Marche, 2014, Figure 20 and significantly better than the results of the Saint-Venant solver (see also e.g. Figure 18 of Hou et al, 2013 and Figure 19 of Nikolos and Delis, 2009) which have too sharp fronts.

We adapt the mesh based on the wavelet-estimated error on the free surface elevation.

```
#if TREE
event adapt (i++) {
scalar η[];
foreach()
η[] = h[] > dry ? h[] + zb[] : 0;
boundary ({η});
astats s = adapt_wavelet ({η}, (double[]){3e-4}, MAXLEVEL, MINLEVEL);
fprintf (stderr, "# refined %d cells, coarsened %d cells\n", s.nf, s.nc);
}
#endif
```