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 2562 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.

Evolution of free surface (left) and mesh (right)

Evolution of free surface (left) and mesh (right)

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.

Timeseries of surface elevation. Numerical results and experimental data (symbols).

Timeseries of surface elevation. Numerical results and experimental data (symbols).

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