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 $h$ so that the water
surface elevation $\eta$ 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)
      foreach_dimension()
	g.x = gradient (zb[-1] + h[-1], eta, zb[1] + h[1])/4.;
    else
      foreach_dimension()
	g.x = (zb[1] - zb[-1])/(2.*Delta);
    // reconstruct water depth h from eta and zb
    foreach_child() {
      double etac = eta;
      foreach_dimension()
	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 $\eta$. We need to do this
    because we cannot assume a priori that the surrounding wet cells are
    necessarily close to e.g. $\eta = 0$. */

    double v = 0., eta = 0.; // water surface elevation
    // 3x3 neighbourhood
    foreach_neighbor(1)
      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 $\eta$ (of the
    parent cell i.e. a first-order interpolation in contrast to the wet
    case above) and $z_b$ of the child cells. */
    
    // reconstruct water depth h from eta and zb
    foreach_child()
      h[] = max(0., eta - zb[]);
  }
}

/**
Cell restriction is simpler. We first compute the depth-weighted
average of $\eta$ over all the children... */

static void restriction_elevation (Point point, scalar h)
{
  double eta = 0., v = 0.;
  foreach_child()
    if (h[] > dry) {
      eta += h[]*(zb[] + h[]);
      v += h[];
    }

  /**
  ... and use this in combination with $z_b$ (of the coarse cell) to
  compute the water depth $h$.  */
    
  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;
  foreach_neighbor(1)
    if (h[] <= dry) {
      wet = false;
      break;
    }
  if (wet)
    refine_linear (point, h);
  else {
    double hc = h[], zc = zb[];
    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;
  h.prolongation = prolongation_elevation;
  h.restriction = restriction_elevation;
  h.dirty = true;
}
#else // Cartesian
void conserve_elevation (void) {}
#endif