# A solver for the Saint-Venant equations

The Saint-Venant equations can be written in integral form as the hyperbolic system of conservation laws ${\partial }_{t}{\int }_{\Omega }\mathbf{\text{q}}d\Omega ={\int }_{\partial \Omega }\mathbf{\text{f}}\left(\mathbf{\text{q}}\right)\cdot \mathbf{\text{n}}d\partial \Omega -{\int }_{\Omega }hg\nabla {z}_{b}$ where $\Omega$ is a given subset of space, $\partial \Omega$ its boundary and $\mathbf{\text{n}}$ the unit normal vector on this boundary. For conservation of mass and momentum in the shallow-water context, $\Omega$ is a subset of bidimensional space and $\mathbf{\text{q}}$ and $\mathbf{\text{f}}$ are written $\mathbf{\text{q}}=\left(\begin{array}{c}\hfill h\hfill \\ \hfill h{u}_{x}\hfill \\ \hfill h{u}_{y}\hfill \end{array}\right),\phantom{\rule{0.278em}{0ex}}\phantom{\rule{0.278em}{0ex}}\phantom{\rule{0.278em}{0ex}}\phantom{\rule{0.278em}{0ex}}\phantom{\rule{0.278em}{0ex}}\phantom{\rule{0.278em}{0ex}}\mathbf{\text{f}}\left(\mathbf{\text{q}}\right)=\left(\begin{array}{cc}\hfill h{u}_{x}\hfill & \hfill h{u}_{y}\hfill \\ \hfill h{u}_{x}^{2}+\frac{1}{2}g{h}^{2}\hfill & \hfill h{u}_{x}{u}_{y}\hfill \\ \hfill h{u}_{x}{u}_{y}\hfill & \hfill h{u}_{y}^{2}+\frac{1}{2}g{h}^{2}\hfill \end{array}\right)$ where $\mathbf{\text{u}}$ is the velocity vector, $h$ the water depth and ${z}_{b}$ 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 ${z}_{b}$ and the flow speed $\mathbf{\text{u}}$. $\eta$ is the water level i.e. ${z}_{b}+h$. Note that the order of the declarations is important as ${z}_{b}$ needs to be refined before $h$ and $h$ before $\eta$.

``````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 $\mathbf{\text{u}}$) are not the conserved variables $h$ and $h\mathbf{\text{u}}$.

``````trace
static void advance_saint_venant (scalar * output, scalar * input,
{
// 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]),
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 $\eta$ is maintained as ${z}_{b}+h$ whenever cells are refined or restrictioned.

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

static void restriction_eta (Point point, scalar eta)
{
η[] = 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 $h\mathbf{\text{u}}$ 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}) {
#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$, $\mathbf{\text{u}}$, ${z}_{b}$. The “interface” topography ${z}_{lr}$ 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});``````

We store the divergence of the fluxes in the update fields. Note that these are updates for $h$ and $h\mathbf{\text{u}}$ (not $\mathbf{\text{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) ${S}_{g}=h\left(\begin{array}{c}\hfill 0\hfill \\ \hfill \frac{g}{2}h{\partial }_{\lambda }{m}_{\theta }+{f}_{G}{u}_{y}\hfill \\ \hfill \frac{g}{2}h{\partial }_{\theta }{m}_{\lambda }-{f}_{G}{u}_{x}\hfill \end{array}\right)$ with ${f}_{G}={u}_{y}{\partial }_{\lambda }{m}_{\theta }-{u}_{x}{\partial }_{\theta }{m}_{\lambda }$

``````      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 = qmalloc (nl, 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"``````