saint-venant-implicit.h

Go to the documentation of this file.

/** @file saint-venant-implicit.h
 */
/**
# A semi-implicit Saint-Venant solver
 
<div class="message">
Note that the [multilayer solver](layered/hydro.h) provides the same
functionality and should be prefered for most applications.</div>
 
We solve the [Saint-Venant equations](saint-venant.h) semi-implicitly
to lift the timestep restriction due to gravity waves. This is just a
particular application of the "all Mach" semi-implicit solver. We will
use the Bell-Collela-Glaz advection scheme to transport mass and
momentum. */
 
#include "all-mach.h"
#include "tracer.h"
 
/**
In addition to the momentum \f$\mathbf{q}=h\mathbf{u}\f$ defined by the
all Mach solver, we will need the fluid depth *h* (i.e. the density
\f$\rho\f$) and the topography *zb*. Both the momentum \f$\mathbf{q}\f$ and
mass \f$h\f$ are advected tracers. The acceleration of gravity is
*G*. *dry* is the fluid depth below which a cell is considered
"dry". */
 
scalar zb[], h[], * tracers = {q,h};
double G = 1.;
double dry = 1e-4;
 
/**
We need fields to store the (varying) state fields \f$\rho c^2\f$,
\f$\alpha\f$ and the acceleration \f$\mathbf{a}\f$. */
 
scalar rhoc2v[];
vector alphav[], av[];
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) {
  rho = h;
  alpha = alphav;
  a = av;
  rhoc2 = rhoc2v;
 
  /**
  We set limiting for *q*, *h*. */
  
  for (int _s = 0; _s < 1; _s++) /* scalar in {q,h} */
    s.gradient = minmod2;
 
  /**
  As well as default values. */
 
  for (int _i = 0; _i < _N; _i++) /* foreach */
    h[] = 1.;
 
  /**
  On trees, we ensure that limiting is also applied to prolongation. */
  
  #if TREE
  for (int _s = 0; _s < 1; _s++) /* scalar in {q,h} */
    s.prolongation = refine_linear;
  #endif
 
  /**
  We setup the default display. */
 
  display ("squares (color = 'h > 0 ? zb + h : nodata', spread = -1);");
}
 
/**
After user initialisation we set the initial pressure and apply
boundary conditions. The boundary conditions for \f$\rho=h\f$ and
\f$\mathbf{q}\f$ are already applied by the [all Mach
solver](all-mach.h). */
 
/** @brief Event: init (i = 0) */
void event_init (void) {
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    p[] = G*sq(h[])/2.;
    /* dim: h[] == Delta */;
  }
}
 
/**
By default the timestep is only limited by the CFL condition on the
advection velocity field. We add the option to define an "acoustic CFL
number" which takes into account the speed of gravity waves. */
 
double CFLa = HUGE;
 
/** @brief Event: stability (i++) */
void event_stability (void) {
  if (CFLa < HUGE)
    foreach (reduction (min:dtmax))
      if (h[] > dry) {
  double dt = CFLa*Delta/sqrt(G*h[]);
  if (dt < dtmax)
    dtmax = dt;
      }
}
 
/**
The properties of the fluid are given by the "Saint-Venant equation of
state" i.e. \f$p = gh^2/2\f$, \f$\rho c^2 = gh^2\f$. */
 
/** @brief Event: properties (i++) */
void event_properties (void) {
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    rhoc2v[] = G*sq(max(h[],dry));
    ps[] = rhoc2v[]/2.;
  }
 
  /**
  The specific volume \f$\alpha=1/\rho\f$ is constructed by averaging,
  taking into account dry states. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
    if ((h[] > dry && h[-1] > dry) ||
  (h[] > dry && h[] + zb[] >= zb[-1]) ||
  (h[-1] > dry && h[-1] + zb[-1] >= zb[]))
      alphav.x[] = 2./(max(h[],dry) + max(h[-1],dry));
    else
      alphav.x[] = 0.;
  }
}
 
/**
## Topographic source term
 
The acceleration due to the topography is \f$- g\nabla z_b\f$. On regular
Cartesian grids we can simply do */
 
#if !TREE
/** @brief Event: acceleration (i++) */
void event_acceleration (void) {
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    if (alpha.x[])
      av.x[] -= G*(zb[] - zb[-1])/Delta;
}
#else // TREE
 
/**
On trees things are a bit more complicated due to the necessity to
verify the "lake-at-rest" condition. For a lake, the pressure gradient
must balance the topographic source term i.e.
\f[
\frac{1}{h}\nabla p = a = -g\nabla z_b
\f]
with \f$a\f$ the acceleration. Using the equation of state and introducing
\f[
\eta = z_b + h
\f]
the elevation of the free surface (constant for a lake at rest), we get
\f[
\frac{1}{2h}\nabla gh^2 = -g\nabla (\eta - h) = g\nabla h
\f]
This identity is obviously verified mathematically, however it is not
necessarily verified by the discrete gradient operator. In the case of
Cartesian meshes it is simple to show that the naive discrete gradient
operator we used above verifies this identity. For tree meshes
this is not generally the case due to the prolongation operator used
to fill ghost cells at refinement boundaries. Rather than trying to
redefine the prolongation operator, we discretise the topographic
source term as
\f[
a = -g\nabla z_b = -g \nabla\eta + \frac{g}{2h}\nabla h^2
\f]
This is strictly identical mathematically, however if we now use the
same discrete operator to estimate \f$\nabla p\f$ and \f$\nabla h^2\f$, the
*discrete* lake-at-rest condition becomes 
\f[ \nabla\eta = 0 \f]
which is much easier to verify.
 
To do so, we define the following restriction and prolongation
operators for \f$\eta\f$. The main idea is to avoid using the elevation of
dry cells. Restriction is done by averaging the elevation of wet
cells. */
 
void eta_restriction (Point point, scalar s)
{
  double sum = 0., n = 0.;
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    if (h[] > dry) {
      sum += s[];
      n++;
    }
  s[] = n > 0 ? sum/n : nodata;
}
 
/**
Prolongation uses linear interpolation with a cell-centered gradient
computed from simple finite-differences of wet cells. */
 
void eta_prolongation (Point point, scalar s)
{
  coord g;
  for (int _d = 0; _d < dimension; _d++) {
    if (h[] > dry) {
      if (h[1] > dry && h[-1] > dry)
  g.x = (s[1] - s[-1])/2.;
      else if (h[1] > dry)
  g.x = s[1] - s[];
      else if (h[-1] > dry)
  g.x = s[] - s[-1];
      else
  g.x = 0.;
    }
    else
      g.x = 0.;
  }
  double sc = s[];
  for (int _c = 0; _c < 4; _c++) /* foreach_child */ {
    s[] = sc;
    for (int _d = 0; _d < dimension; _d++)
      s[] += child.x*g.x/4.;
  }
}
 
/**
One can check that the prolongation values constructed using these
functions verify \f$\eta = constant\f$ for a lake-at-rest. */
 
/** @brief Event: acceleration (i++) */
void event_acceleration (void) {
  
  /**
  To compute the acceleration due to topography, we first compute
  \f$\eta\f$ and \f$h^2\f$. */
  
  scalar eta[], h2[];
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    h2[] = sq(h[]);
    eta[] = zb[] + h[];
  }
 
  /**
  We then make sure that \f$\eta\f$ uses our restriction/prolongation
  functions. \f$h^2\f$ uses the same prolongation functions as \f$p\f$ by
  default. */
 
  eta.prolongation = eta_prolongation;
  eta.restriction = eta_restriction;
    
  /**
  We then compute the acceleration as
  \f[
  a = -g \nabla\eta + g \frac{\alpha}{2}\nabla h^2  
  \f]
  */
  
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    if (alpha.x[])
      av.x[] += G*(eta[-1] - eta[] + alpha.x[]/2.*(h2[] - h2[-1]))/Delta;
}
 
#endif // TREE
 
/**
The utility functions in *elevation.h* need to know which gradient we
are using. */
 
double (* gradient)  (double, double, double) = minmod2;
 
#include "elevation.h"
#include "gauges.h"