elevation.h

Go to the documentation of this file.

/** @file elevation.h
 */
/**
# Conservation of water surface elevation 
 
When using the default adaptive reconstruction of variables, the
[Saint-Venant solver](saint-venant.h) or the [layered solver](hydro.h)
will conserve the water depth when cells are refined or
coarsened. However, this will not necessarily ensure that the
"lake-at-rest" condition (i.e. a constant water surface elevation) is
also preserved. In what follows, we redefine the *prolongation()* and
*restriction()* methods of the water depth \f$h\f$ so that the water
surface elevation \f$\eta\f$ is conserved.
 
We start with the reconstruction of fine "wet" cells: */
 
#if TREE
static double default_sea_level = 0.;
 
static void refine_elevation (Point point, scalar h)
{
  // reconstruction of fine cells using elevation (rather than water depth)
  // (default refinement conserves mass but not lake-at-rest)
  if (h[] >= dry) {
    double eta = zb[] + h[];   // water surface elevation  
    coord g; // gradient of eta
    if (gradient)
      for (int _d = 0; _d < dimension; _d++)
  g.x = gradient (zb[-1] + h[-1], eta, zb[1] + h[1])/4.;
    else
      for (int _d = 0; _d < dimension; _d++)
  g.x = (zb[1] - zb[-1])/(2.*Delta);
    // reconstruct water depth h from eta and zb
    for (int _c = 0; _c < 4; _c++) /* foreach_child */ {
      double etac = eta;
      for (int _d = 0; _d < dimension; _d++)
  etac += g.x*child.x;
      h[] = max(0., etac - zb[]);
    }
  }
  else {
 
    /**
    The "dry" case is a bit more complicated. We look in a 3x3
    neighborhood of the coarse parent cell and compute a depth-weighted
    average of the "wet" surface elevation \f$\eta\f$. We need to do this
    because we cannot assume a priori that the surrounding wet cells are
    necessarily close to e.g. \f$\eta = 0\f$. */
 
    double v = 0., eta = 0.; // water surface elevation
    // 3x3 neighbourhood
    for (int _n = 0; _n < 1; _n++)
      if (h[] >= dry) {
  eta += h[]*(zb[] + h[]);
  v += h[];
      }
    if (v > 0.)
      eta /= v; // volume-averaged eta of neighbouring wet cells
    else
 
      /**
      If none of the surrounding cells is wet, we set a default sealevel. */
 
      eta = default_sea_level;
 
    /**
    We then reconstruct the water depth in each child using \f$\eta\f$ (of the
    parent cell i.e. a first-order interpolation in contrast to the wet
    case above) and \f$z_b\f$ of the child cells. */
    
    // reconstruct water depth h from eta and zb
    for (int _c = 0; _c < 4; _c++) /* foreach_child */
      h[] = max(0., eta - zb[]);
  }
}
 
/**
Cell restriction is simpler. We first compute the depth-weighted
average of \f$\eta\f$ over all the children... */
 
static void restriction_elevation (Point point, scalar h)
{
  double eta = 0., v = 0.;
  for (int _c = 0; _c < 4; _c++) /* foreach_child */
    if (h[] > dry) {
      eta += h[]*(zb[] + h[]);
      v += h[];
    }
 
  /**
  ... and use this in combination with \f$z_b\f$ (of the coarse cell) to
  compute the water depth \f$h\f$.  */
    
  if (v > 0.)
    h[] = max(0., eta/v - zb[]);
  else // dry cell
    h[] = 0.;
}
 
/**
We also need to define a consistent prolongation function. For cells
which are entirely surrounded by wet cells, we can use the standard
linear refinement function, otherwise we use straight injection from
the parent cell. */
 
static void prolongation_elevation (Point point, scalar h)
{
  bool wet = true;
  for (int _n = 0; _n < 1; _n++)
    if (h[] <= dry) {
      wet = false;
      break;
    }
  if (wet)
    refine_linear (point, h);
  else {
    double hc = h[], zc = zb[];
    for (int _c = 0; _c < 4; _c++) /* foreach_child */ {
      h[] = hc;
      zb[] = zc;
    }
  }
}
 
/**
Finally we define a function which will be called by the user to apply
these reconstructions.  */
 
void conserve_elevation (void)
{
  h.refine  = refine_elevation;
  set_prolongation (h, prolongation_elevation);
  set_restriction (h, restriction_elevation);
}
#else // Cartesian
void conserve_elevation (void) {}
#endif