src/green-naghdi.h

A solver for the Green-Naghdi equations

The Green-Naghdi equations (also known as the Serre or fully nonlinear Boussinesq equations) can be seen as an extension of the Saint-Venant equations to the next order O(μ2) in term of the shallowness parameter μ=h02L2 with h0 the typical depth and L the typical horizontal scale. In contrast to the Saint-Venant equations the Green-Naghdi equations have dispersive wave solutions.

A more detailed description of the context and numerical scheme implemented here is given in Popinet, 2015.

The solver is built by adding a source term to the momentum equation of the Saint-Venant solver. Following Bonneton et al, 2011, this source term can be written t(hu)+=h(gαdηD) where D verifies αdh𝒯(D)+hD=b and b=[gαdη+𝒬1(u)] With 𝒯 a linear operator to be defined below, as well as 𝒬1(u).

Before including the Saint-Venant solver, we need to overload the default update function of the predictor-corrector scheme in order to add our source term.

#include "predictor-corrector.h"

static double update_green_naghdi (scalar * current, scalar * updates,
				   double dtmax);

event defaults (i = 0)
  update = update_green_naghdi;

#include "saint-venant.h"

The linear system can be inverted with the multigrid Poisson solver. We declare D as a global variable so that it can be re-used as initial guess for the Poisson solution. The solver statistics will be stored in mgD. The breaking parameter defines the slope above which dispersive terms are turned off. The αd parameter controls the optimisation of the dispersion relation (see Bonneton et al, 2011).

#include "poisson.h"

vector D[];
mgstats mgD;
double breaking = 1., alpha_d = 1.153;

We first define some useful macros, following the notations in Bonneton et al, 2011.

Simple centered-difference approximations of the first and second derivatives of a field.

#define dx(s)  ((s[1,0] - s[-1,0])/(2.*Δ))
#define dy(s)  ((s[0,1] - s[0,-1])/(2.*Δ))
#define d2x(s) ((s[1,0] + s[-1,0] - 2.*s[])/sq(Δ))
#define d2y(s) ((s[0,1] + s[0,-1] - 2.*s[])/sq(Δ))
#define d2xy(s) ((s[1,1] - s[1,-1] - s[-1,1] + s[-1,-1])/sq(2.*Δ))

The definitions of the 1 and 2 operators 1[h,zb]w=13h(h3w)h2wzb=h[h3w+w(h+12zb)]2[h,zb]w=12h(h2w)+wzb=h2w+w(zb+h)

#define R1(h,zb,w) (-h[]*(h[]/3.*dx(w) + w[]*(dx(h) + dx(zb)/2.)))
#define R2(h,zb,w) (h[]/2.*dx(w) + w[]*(dx(zb) + dx(h)))

Inversion of the linear system

To invert the linear system which defines D, we need to write discretised versions of the residual and relaxation operators. The functions take generic multilevel parameters and a user-defined structure which contains solution-specific parameters, in our case a list of the h, zb and wet fields.

static double residual_GN (scalar * a, scalar * r, scalar * resl, void * data)
{

We first recover all the parameters from the generic pointers and rename them according to our notations.

  scalar * list = data;
  scalar h = list[0], zb = list[1], wet = list[2];
  vector D = vector(a[0]), b = vector(r[0]), res = vector(resl[0]);

The general form for 𝒯 is 𝒯[h,zb]w=1[h,zb](w)+2[h,zb](zbw) which gives the linear problem αdh𝒯(D)+hD=b αdh1(D)+αdh2(zbD)+hD=b αdh2[h3(D)+(D)(h+12zb)]+αdh[h2(zbD)+(zbD)η]+hD=b Expanding the operators for the x-component gives αd3x(h3xDx)+h[αd(xηxzb+h2x2zb)+1]Dx+αdh[(h2xy2zb+xηyzb)Dy+h2yzbxDyh23xy2DyhyDy(xh+12xzb)]=bx The y-component is obtained by rotation of the indices.

  double maxres = 0.;
  foreach (reduction(max:maxres))
    foreach_dimension() {
      if (wet[-1] == 1 && wet[] == 1 && wet[1] == 1) {
	double hc = h[], dxh = dx(h), dxzb = dx(zb), dxeta = dxzb + dxh;
	double hl3 = (hc + h[-1])/2.; hl3 = cube(hl3);
	double hr3 = (hc + h[1])/2.;  hr3 = cube(hr3);

Finally we translate the formula above to its discrete version.

	res.x[] = b.x[] -
	  (-alpha_d/3.*(hr3*D.x[1] + hl3*D.x[-1] - 
			(hr3 + hl3)*D.x[])/sq(Δ) +
	   hc*(alpha_d*(dxeta*dxzb + hc/2.*d2x(zb)) + 1.)*D.x[] +
	   alpha_d*hc*((hc/2.*d2xy(zb) + dxeta*dy(zb))*D.y[] + 
		       hc/2.*dy(zb)*dx(D.y) - sq(hc)/3.*d2xy(D.y)
		       - hc*dy(D.y)*(dxh + dxzb/2.)));

The function also need to return the maximum residual.

	if (fabs (res.x[]) > maxres)
	  maxres = fabs (res.x[]);
      }
      else
	res.x[] = 0.;
    }
  boundary (resl);

The maximum residual is normalised by gravity i.e. the tolerance is the relative acceleration times the depth.

  return maxres/G;
}

The relaxation function is built by copying and pasting the residual implementation above and inverting manually for Dx.

static void relax_GN (scalar * a, scalar * r, int l, void * data)
{
  scalar * list = data;
  scalar h = list[0], zb = list[1], wet = list[2];
  vector D = vector(a[0]), b = vector(r[0]);
  foreach_level_or_leaf (l)
    foreach_dimension() {
      if (h[] > dry && wet[-1] == 1 && wet[] == 1 && wet[1] == 1) {
	double hc = h[], dxh = dx(h), dxzb = dx(zb), dxeta = dxzb + dxh;
	double hl3 = (hc + h[-1])/2.; hl3 = cube(hl3);
	double hr3 = (hc + h[1])/2.;  hr3 = cube(hr3);
	D.x[] = (b.x[] -
		 (-alpha_d/3.*(hr3*D.x[1] + hl3*D.x[-1])/sq(Δ) +
		  alpha_d*hc*((hc/2.*d2xy(zb) + dxeta*dy(zb))*D.y[] + 
			      hc/2.*dy(zb)*dx(D.y) - sq(hc)/3.*d2xy(D.y)
			      - hc*dy(D.y)*(dxh + dxzb/2.))))/
	  (alpha_d*(hr3 + hl3)/(3.*sq(Δ)) + 
	   hc*(alpha_d*(dxeta*dxzb + hc/2.*d2x(zb)) + 1.));
      }
      else
	D.x[] = 0.;
    }
}

Source term computation

To add the source term to the Saint-Venant system we overload the default update function with this one. The function takes a list of the current evolving scalar fields and fills the corresponding updates with the source terms.

static double update_green_naghdi (scalar * current, scalar * updates,
				   double dtmax)
{
  double dt = update_saint_venant (current, updates, dtmax);
  scalar h = current[0];
  vector u = vector(current[1]);

We first compute the right-hand-side b. The general form for the 𝒬1 operator is (eq. (9) of Bonneton et al, 2011). 𝒬1[h,zb](u)=21(c)+2(d) with c=1u2u+(u)2=xuxyuy+xuyyux+(xux+yuy)2d=u(u)zb=ux(uxx2zb+uyxy2zb)+uy(uyy2zb+uxxy2zb)=ux2x2zb+uy2y2zb+2uxuyxy2zb Note that we declare field c and d in a new scope, so that the corresponding memory is freed after we have computed b.

  vector b[];
  {
    scalar c[], d[];
    foreach() {
      double dxux = dx(u.x), dyuy = dy(u.y);
      c[] = - dxux*dyuy + dx(u.y)*dy(u.x) + sq(dxux + dyuy);
      d[] = sq(u.x[])*d2x(zb) + sq(u.y[])*d2y(zb) + 2.*u.x[]*u.y[]*d2xy(zb);
    }
    boundary ({c,d});

We can now compute b b=[gαdη+𝒬1(u)]

    foreach()
      foreach_dimension()
        b.x[] = h[]*(G/alpha_d*dx(η) - 2.*R1(h,zb,c) + R2(h,zb,d));
  }

We declare a new field which will track whether cells are completely wet. This is necessary to turn off dispersive terms close to the shore so that lake-at-rest balance is maintained.

  scalar wet[];
  foreach()
    wet[] = h[] > dry;
  boundary ({wet});

Finally we solve the linear problem with the multigrid solver using the relaxation and residual functions defined above. We need to restrict h and zb as their values will be required for relaxation on coarse grids.

  scalar * list = {h, zb, wet};
  restriction (list);
  mgD = mg_solve ((scalar *){D}, (scalar *){b},
		  residual_GN, relax_GN, list, mgD.nrelax);

We can then compute the updates for hu.

  vector dhu = vector(updates[1]);

We only apply the Green-Naghdi source term when the slope of the free surface is smaller than breaking. The idea is to turn off dispersion in areas where the wave is “breaking” (i.e. in hydraulic jumps). We also turn off dispersive terms close to shore so that lake-at-rest balance is maintained.

  foreach()
    if (fabs(dx(η)) < breaking && fabs(dy(η)) < breaking)
      foreach_dimension()
	if (wet[-1] == 1 && wet[] == 1 && wet[1] == 1)
	  dhu.x[] += h[]*(G/alpha_d*dx(η) - D.x[]);

  return dt;
}

Usage

Examples

Tests