src/saint-venant.h

A solver for the Saint-Venant equations

The Saint-Venant equations can be written in integral form as the hyperbolic system of conservation laws tΩqdΩ=Ωf(q)ndΩΩhgzb where Ω is a given subset of space, Ω its boundary and n the unit normal vector on this boundary. For conservation of mass and momentum in the shallow-water context, Ω is a subset of bidimensional space and q and f are written q=(hhuxhuy),f(q)=(huxhuyhux2+12gh2huxuyhuxuyhuy2+12gh2) where u is the velocity vector, h the water depth and zb the height of the topography. See also Popinet, 2011 for a more detailed introduction.

User variables and parameters

The primary fields are the water depth h, the bathymetry zb and the flow speed u. η is the water level i.e. zb+h. Note that the order of the declarations is important as zb needs to be refined before h and h before η.

scalar zb[], h[], η[];
vector u[];

The only physical parameter is the acceleration of gravity G. Cells are considered “dry” when the water depth is less than the dry parameter (this should not require tweaking).

double G = 1.;
double dry = 1e-10;

By default there is only a single layer i.e. this is the classical Saint-Venant system above. This can be changed by setting nl to a different value. The extension of the Saint-Venant system to multiple layers is implemented in multilayer.h.

int nl = 1;
#include "multilayer.h"

Time-integration

Setup

Time integration will be done with a generic predictor-corrector scheme.

#include "predictor-corrector.h"

The generic time-integration scheme in predictor-corrector.h needs to know which fields are updated. The list will be constructed in the defaults event below.

scalar * evolving = NULL;

We need to overload the default advance function of the predictor-corrector scheme, because the evolving variables (h and u) are not the conserved variables h and hu.

trace
static void advance_saint_venant (scalar * output, scalar * input, 
				  scalar * updates, double dt)
{
  // recover scalar and vector fields from lists
  scalar hi = input[0], ho = output[0], dh = updates[0];
  vector * uol = (vector *) &output[1];

  // new fields in ho[], uo[]
  foreach() {
    double hold = hi[];
    ho[] = hold + dt*dh[];
    η[] = zb[] + ho[];
    if (ho[] > dry) {
      for (int l = 0; l < nl; l++) {
        vector uo = vector(output[1 + dimension*l]);
      	vector ui = vector(input[1 + dimension*l]),
	  dhu = vector(updates[1 + dimension*l]);
	foreach_dimension()
	  uo.x[] = (hold*ui.x[] + dt*dhu.x[])/ho[];
      }

In the case of multiple layers we add the viscous friction between layers.

      if (nl > 1)
	vertical_viscosity (point, ho[], uol, dt);
    }
    else // dry
      for (int l = 0; l < nl; l++) {
        vector uo = vector(output[1 + dimension*l]);
        foreach_dimension()
	  uo.x[] = 0.;
      }
  }
    
  // fixme: on trees eta is defined as eta = zb + h and not zb +
  // ho in the refine_eta() and restriction_eta() functions below
  scalar * list = list_concat ({ho, η}, (scalar *) uol);
  boundary (list);
  free (list);
}

When using an adaptive discretisation (i.e. a tree)., we need to make sure that η is maintained as zb+h whenever cells are refined or restrictioned.

#if TREE
static void refine_eta (Point point, scalar η)
{
  foreach_child()
    η[] = zb[] + h[];
}

static void restriction_eta (Point point, scalar η)
{
  η[] = zb[] + h[];
}
#endif

Computing fluxes

Various approximate Riemann solvers are defined in riemann.h.

#include "riemann.h"

trace
double update_saint_venant (scalar * evolving, scalar * updates, double dtmax)
{

We first recover the currently evolving height and velocity (as set by the predictor-corrector scheme).

  scalar h = evolving[0], dh = updates[0];
  vector u = vector(evolving[1]);

Fh and Fq will contain the fluxes for h and hu respectively and S is necessary to store the asymmetric topographic source term.

  face vector Fh[], S[];
  tensor Fq[];

The gradients are stored in locally-allocated fields. First-order reconstruction is used for the gradient fields.

  vector gh[], geta[];
  tensor gu[];
  for (scalar s in {gh, geta, gu}) {
    s.gradient = zero;
    #if TREE
      s.prolongation = refine_linear;
    #endif
  }
  gradients ({h, η, u}, {gh, geta, gu});

We go through each layer.

  for (int l = 0; l < nl; l++) {

We recover the velocity field for the current layer and compute its gradient (for the first layer the gradient has already been computed above).

    vector u = vector (evolving[1 + dimension*l]);
    if (l > 0)
      gradients ((scalar *) {u}, (vector *) {gu});

The faces which are “wet” on at least one side are traversed.

    foreach_face (reduction (min:dtmax)) {
      double hi = h[], hn = h[-1];
      if (hi > dry || hn > dry) {

Left/right state reconstruction

The gradients computed above are used to reconstruct the left and right states of the primary fields h, u, zb. The “interface” topography zlr is reconstructed using the hydrostatic reconstruction of Audusse et al, 2004

	double dx = Δ/2.;
	double zi = η[] - hi;
	double zl = zi - dx*(geta.x[] - gh.x[]);
	double zn = η[-1] - hn;
	double zr = zn + dx*(geta.x[-1] - gh.x[-1]);
	double zlr = max(zl, zr);
	
	double hl = hi - dx*gh.x[];
	double up = u.x[] - dx*gu.x.x[];
	double hp = max(0., hl + zl - zlr);
	
	double hr = hn + dx*gh.x[-1];
	double um = u.x[-1] + dx*gu.x.x[-1];
	double hm = max(0., hr + zr - zlr);

Riemann solver

We can now call one of the approximate Riemann solvers to get the fluxes.

	double fh, fu, fv;
	kurganov (hm, hp, um, up, Δ*cm[]/fm.x[], &fh, &fu, &dtmax);
	fv = (fh > 0. ? u.y[-1] + dx*gu.y.x[-1] : u.y[] - dx*gu.y.x[])*fh;

Topographic source term

In the case of adaptive refinement, care must be taken to ensure well-balancing at coarse/fine faces (see notes/balanced.tm).

        #if TREE
	if (is_prolongation(cell)) {
	  hi = coarse(h);
	  zi = coarse(zb);
	}
	if (is_prolongation(neighbor(-1))) {
	  hn = coarse(h,-1);
	  zn = coarse(zb,-1);
	}
        #endif
	
	double sl = G/2.*(sq(hp) - sq(hl) + (hl + hi)*(zi - zl));
	double sr = G/2.*(sq(hm) - sq(hr) + (hr + hn)*(zn - zr));

Flux update

	Fh.x[]   = fm.x[]*fh;
	Fq.x.x[] = fm.x[]*(fu - sl);
	S.x[]    = fm.x[]*(fu - sr);
	Fq.y.x[] = fm.x[]*fv;
      }
      else // dry
	Fh.x[] = Fq.x.x[] = S.x[] = Fq.y.x[] = 0.;
    }

    boundary_flux ({Fh, S, Fq});

Updates for evolving quantities

We store the divergence of the fluxes in the update fields. Note that these are updates for h and hu (not u).

    vector dhu = vector(updates[1 + dimension*l]);
    foreach() {
      double dhl =
	layer[l]*(Fh.x[1,0] - Fh.x[] + Fh.y[0,1] - Fh.y[])/(cm[]*Δ);
      dh[] = - dhl + (l > 0 ? dh[] : 0.);
      foreach_dimension()
	dhu.x[] = (Fq.x.x[] + Fq.x.y[] - S.x[1,0] - Fq.x.y[0,1])/(cm[]*Δ);

For multiple layers we need to store the divergence in each layer.

      if (l < nl - 1) {
	scalar div = wl[l];
	div[] = dhl;
      }

We also need to add the metric terms. They can be written (see eq. (8) of Popinet, 2011) Sg=h(0g2hλmθ+fGuyg2hθmλfGux) with fG=uyλmθuxθmλ

      double dmdl = (fm.x[1,0] - fm.x[])/(cm[]*Δ);
      double dmdt = (fm.y[0,1] - fm.y[])/(cm[]*Δ);
      double fG = u.y[]*dmdl - u.x[]*dmdt;
      dhu.x[] += h[]*(G*h[]/2.*dmdl + fG*u.y[]);
      dhu.y[] += h[]*(G*h[]/2.*dmdt - fG*u.x[]);
    }
  }

For multiple layers we need to add fluxes between layers.

  if (nl > 1)
    vertical_fluxes ((vector *) &evolving[1], (vector *) &updates[1], wl, dh);
    
  return dtmax;
}

Initialisation and cleanup

We use the main time loop (in the predictor-corrector scheme) to setup the initial defaults.

event defaults (i = 0)
{
  assert (ul == NULL && wl == NULL);
  assert (nl > 0);
  ul = vectors_append (ul, u);
  for (int l = 1; l < nl; l++) {
    scalar w = new scalar;
    vector u = new vector;
    ul = vectors_append (ul, u);
    wl = list_append (wl, w);
  }

  evolving = list_concat ({h}, (scalar *) ul);
  foreach()
    for (scalar s in evolving)
      s[] = 0.;
  boundary (evolving);

By default, all the layers have the same relative thickness.

  layer = malloc (nl*sizeof(double));
  for (int l = 0; l < nl; l++)
    layer[l] = 1./nl;

We overload the default ‘advance’ and ‘update’ functions of the predictor-corrector scheme and setup the prolongation and restriction methods on trees.

  advance = advance_saint_venant;
  update = update_saint_venant;
#if TREE
  for (scalar s in {h,zb,u,η}) {
    s.refine = s.prolongation = refine_linear;
    s.restriction = restriction_volume_average;
  }
  η.refine  = refine_eta;
  η.restriction = restriction_eta;
#endif
}

The event below will happen after all the other initial events to take into account user-defined field initialisations.

event init (i = 0)
{
  foreach()
    η[] = zb[] + h[];
  boundary (all);
}

At the end of the simulation, we free the memory allocated in the defaults event.

event cleanup (i = end, last) {
  free (evolving);
  free (layer);
  free (ul), ul = NULL;
  free (wl), wl = NULL;
}

#include "elevation.h"

Usage

Examples

Tests