sandbox/M1EMN/Exemples/saint-venantB.h
A solver for the Saint-Venant equations
This is the original file, only one change is done (multilayer-multilayerB) The Saint-Venant equations can be written in integral form as the hyperbolic system of conservation laws \displaystyle \partial_t \int_{\Omega} \mathbf{q} d \Omega = \int_{\partial \Omega} \mathbf{f} ( \mathbf{q}) \cdot \mathbf{n}d \partial \Omega - \int_{\Omega} hg \nabla z_b where \Omega is a given subset of space, \partial \Omega its boundary and \mathbf{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{q} and \mathbf{f} are written \displaystyle \mathbf{q} = \left(\begin{array}{c} h\\ hu_x\\ hu_y \end{array}\right), \;\;\;\;\;\; \mathbf{f} (\mathbf{q}) = \left(\begin{array}{cc} hu_x & hu_y\\ hu_x^2 + \frac{1}{2} gh^2 & hu_xu_y\\ hu_xu_y & hu_y^2 + \frac{1}{2} gh^2 \end{array}\right) where \mathbf{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{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[], eta[];
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 “multilayerB.h”
Multilayer Saint-Venant system with mass exchanges
This is a modification of the original file
The Saint-Venant system is extended to multiple layers following Audusse et al, 2011 as \displaystyle \partial_th + \partial_x\sum_{l=0}^{nl-1}h_lu_l = 0 with \displaystyle h_l = \mathrm{layer}_lh with \mathrm{layer}_l the relative thickness of the layers satisfying \displaystyle \mathrm{layer}_l >= 0,\;\sum_{l=0}^{nl - 1}\mathrm{layer}_l = 1. The momentum equation in each layer is thus \displaystyle \partial_t(h\mathbf{u}_l) + \nabla\cdot\left(h\mathbf{u}_l\otimes\mathbf{u}_l + \frac{gh^2}{2}\mathbf{I}\right) = - gh\nabla z_b + \frac{1}{\mathrm{layer}_l}\left[\mathbf{u}_{l+1/2}G_{l+1/2} - \mathbf{u}_{l-1/2}G_{l-1/2} + \nu\left(\frac{u_{l+1} - u_l}{h_{l+1/2}} - \frac{u_{l} - u_{l-1}}{h_{l-1/2}}\right)\right] where G_{l+1/2} is the relative vertical transport velocity between layers and the second term corresponds to viscous friction between layers.
These last two terms are the only difference with the one layer system.
The horizontal velocity in each layer is stored in ul and the vertical velocity between layers in wl.
vector * ul = NULL;
scalar * wl = NULL;
double * layer;Viscous friction between layers
Boundary conditions on the top and bottom layers need to be added to close the system for the viscous stresses. We chose to impose a Neumann condition on the top boundary i.e. \displaystyle \partial_z u |_t = \dot{u}_t and a Navier slip condition on the bottom i.e. \displaystyle u|_b = u_b + \lambda_b \partial_z u|_b By default the viscosity is zero and we impose free-slip on the top boundary and no-slip on the bottom boundary i.e. \dot{u}_t = 0, \lambda_b = 0, u_b = 0.
\nu B is the yeld stress
double B ;
double nu ;
(const) scalar lambda_b = zeroc, dut = zeroc, u_b = zeroc;For stability, we discretise the viscous friction term implicitly as \displaystyle \frac{(hu_l)_{n + 1} - (hu_l)_{\star}}{\Delta t} = \frac{\nu}{\mathrm{layer}_l} \left( \frac{u_{l + 1} - u_l}{h_{l + 1 / 2}} - \frac{u_l - u_{l - 1}}{h_{l - 1 / 2}} \right)_{n + 1} which can be expressed as the linear system \displaystyle \mathbf{Mu}_{n + 1} = \mathrm{rhs} where \mathbf{M} is a tridiagonal matrix. The lower, principal and upper diagonals are a, b and c respectively.
double nu_eq(double,double);
void vertical_viscosity (Point point, double h, vector * ul, double dt)
{
if (nu == 0.)
return;
double nueq[nl];
double a[nl], b[nl], c[nl], rhs[nl];
foreach_dimension() {The rhs of the tridiagonal system is h_lu_l = h\mathrm{layer}_lu_l.
int l = 0;
for (vector u in ul)
rhs[l] = h*layer[l]*u.x[], l++;*definition of viscosity (it is a function of shear: Bingham flow, and the viscosity depends on pressure in the case of \mu(I):
for example, in the Bingham case, you define the function \displaystyle \tau = \nu (\frac{\partial u}{\partial z} + B) where \nu is the viscosity and \nu B is the yeld stress, we define an equivalent viscosity: \displaystyle \nu_{eq}=\nu (1 + \frac{B}{|\frac{\partial u}{\partial z}|}) Note the regularization introduced to avoid division by zero*
double press = G*h;//eta[]-zb[];
for (l = 1; l < nl; l++) {
vector um = ul[l-1] ;
vector up = ul[l] ;
press -= h*layer[l];
double shear = (h>0? (up.x[]-um.x[])/(h*layer[l]) : HUGE);
nueq[l] = nu_eq(shear,press);
}
nueq[0] = nueq[1];The lower, principal and upper diagonals a, b and c are given by \displaystyle a_{l > 0} = - \left( \frac{\nu \Delta t}{h_{l - 1 / 2}} \right)_{n + 1} \displaystyle c_{l < \mathrm{nl} - 1} = - \left( \frac{\nu \Delta t}{h_{l + 1 / 2}} \right)_{n + 1} \displaystyle b_{0 < l < \mathrm{nl} - 1} = \mathrm{layer}_l h_{n + 1} - a_l - c_l
for (l = 1; l < nl - 1; l++) {
a[l] = - (nueq[l]+nueq[l-1])*dt/(h*(layer[l-1] + layer[l]));
c[l] = - (nueq[l+1]+nueq[l])*dt/(h*(layer[l] + layer[l+1]));
b[l] = layer[l]*h - a[l] - c[l];
}For the top layer the boundary conditions give the (ghost) boundary value \displaystyle u_{\mathrm{nl}} = u_{\mathrm{nl} - 1} + \dot{u}_t h_{\mathrm{nl} - 1}, which gives the diagonal coefficient and right-hand-side \displaystyle b_{\mathrm{nl} - 1} = \mathrm{layer}_{\mathrm{nl} - 1} h_{n + 1} - a_{\mathrm{nl} - 1} \displaystyle \mathrm{rhs}_{\mathrm{nl} - 1} = \mathrm{layer}_{\mathrm{nl} - 1} (hu_{\mathrm{nl} - 1})_{\star} + \nu \Delta t \dot{u}_t
a[nl-1] = - (nueq[l-1]+nueq[l-2])*dt/(h*(layer[nl-2] + layer[nl-1]));
b[nl-1] = layer[nl-1]*h - a[nl-1];
rhs[nl-1] += nueq[l-1]*dt*dut[];For the bottom layer, the boundary conditions give the (ghost) boundary value u_{- 1} \displaystyle u_{- 1} = \frac{2 h_0}{2 \lambda_b + h_0} u_b + \frac{2 \lambda_b - h_0}{2 \lambda_b + h_0} u_0, which gives the diagonal coefficient and right-hand-side \displaystyle b_0 = \mathrm{layer}_0 h_{n + 1} - c_0 + \frac{2 \nu \Delta t}{2 \lambda_b + h_0} \displaystyle \mathrm{rhs}_0 = \mathrm{layer}_0 (hu_0)_{\star} + \frac{2 \nu \Delta t}{2 \lambda_b + h_0} u_b
c[0] = - 2.*dt*nueq[0]/(h*(layer[0] + layer[1]));
b[0] = layer[0]*h - c[0] + 2.*nueq[0]*dt/(2.*lambda_b[] + h*layer[0]);
rhs[0] += 2.*nueq[0]*dt/(2.*lambda_b[] + h*layer[0])*u_b[];We can now solve the tridiagonal system using the Thomas algorithm.
for (l = 1; l < nl; l++) {
b[l] -= a[l]*c[l-1]/b[l-1];
rhs[l] -= a[l]*rhs[l-1]/b[l-1];
}
vector u = ul[nl-1];
u.x[] = a[nl-1] = rhs[nl-1]/b[nl-1];
for (l = nl - 2; l >= 0; l--) {
u = ul[l];
u.x[] = a[l] = (rhs[l] - c[l]*a[l+1])/b[l];
}
}
}Fluxes between layers
The relative vertical velocity between layers l and l+1 is defined as (eq. (2.22) of Audusse et al, 2011) \displaystyle G_{l+1/2} = \sum_{j=0}^{l}(\mathrm{div}_j + \mathrm{layer}_j\mathrm{dh}) with \displaystyle \mathrm{div}_l = \nabla\cdot(h_l\mathbf{u}_l) \displaystyle \mathrm{dh} = - \sum_{l=0}^{nl-1} \mathrm{div}_l
void vertical_fluxes (vector * evolving, vector * updates,
scalar * divl, scalar dh)
{
foreach() {
double Gi = 0., sumjl = 0.;
for (int l = 0; l < nl - 1; l++) {
scalar div = divl[l];
Gi += div[] + layer[l]*dh[];
sumjl += layer[l];
scalar w = div;
w[] = dh[]*sumjl - Gi;
foreach_dimension() {To compute the vertical advection term, we need an estimate of the velocity at l+1/2. This is obtained using simple upwinding according to the sign of the interface velocity \mathrm{Gi} = G_{l+1/2} and the values of the velocity in the l and l+1 layers. Note that the inequality of upwinding is consistent with equs. (5.110) of Audusse et al, 2011 and (77) of Audusse et al, 2011b but not with eq. (2.23) of Audusse et al, 2011.
scalar ub = evolving[l].x, ut = evolving[l + 1].x;
double ui = Gi < 0. ? ub[] : ut[];The flux at l+1/2 is then added to the updates of the bottom layer and substracted from the updates of the top layer.
double flux = Gi*ui;
scalar du_b = updates[l].x, du_t = updates[l + 1].x;
du_b[] += flux/layer[l];
du_t[] -= flux/layer[l + 1];To compute the vertical velocity we use the definition of the mass flux term (eq. 2.13 of Audusse et al, 2011): \displaystyle \mathrm{w}(\mathbf{x},z_{l+1/2}) = \partial_t z_{l+1/2} - G_{l+1/2} + \mathbf{u}_{l+1/2} \cdot \nabla z_{l+1/2} We can write the vertical position of the interface as: \displaystyle z_{l+1/2} = z_{b} + \sum_{j=0}^{l} h_{j} so that the vertical velocity is: \displaystyle \mathrm{w}(\mathbf{x},z_{l+1/2}) = \mathrm{dh}\sum_{j=0}^{l}\mathrm{layer}_{j} - G_{l+1/2} + \mathbf{u}_{l+1/2} \cdot \left[\nabla z_{b} + \nabla h \sum_{j=0}^{l}\mathrm{layer}_{j}\right]
w[] += ui*((zb[1] - zb[-1]) + (h[1] - h[-1])*sumjl)/(2.*Delta);
}
}
}
}——-end of —#include “multilayerB.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{u}) are not the conserved variables h and h\mathbf{u}.
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[];
eta[] = 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, eta}, (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()
eta[] = zb[] + h[];
}
static void restriction_eta (Point point, scalar eta)
{
eta[] = zb[] + h[];
}
#endifComputing 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{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}) {
s.gradient = zero;
#if TREE
s.prolongation = refine_linear;
#endif
}
gradients ({h, eta, 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{u}, z_b. The “interface” topography z_{lr} is reconstructed using the hydrostatic reconstruction of Audusse et al, 2004
double dx = Delta/2.;
double zi = eta[] - hi;
double zl = zi - dx*(geta.x[] - gh.x[]);
double zn = eta[-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, Delta*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 h\mathbf{u} (not \mathbf{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[]*Delta);
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[]*Delta);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) \displaystyle S_g = h \left(\begin{array}{c} 0\\ \frac{g}{2} h \partial_{\lambda} m_{\theta} + f_G u_y\\ \frac{g}{2} h \partial_{\theta} m_{\lambda} - f_G u_x \end{array}\right) with \displaystyle f_G = u_y \partial_{\lambda} m_{\theta} - u_x \partial_{\theta} m_{\lambda}
double dmdl = (fm.x[1,0] - fm.x[])/(cm[]*Delta);
double dmdt = (fm.y[0,1] - fm.y[])/(cm[]*Delta);
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,eta}) {
s.refine = s.prolongation = refine_linear;
s.restriction = restriction_volume_average;
}
eta.refine = refine_eta;
eta.restriction = restriction_eta;
#endif
}The event below will happen after all the other initial events to take into account user-defined field initialisations.
At the end of the simulation, we free the memory allocated in the defaults event.
