src/elevation.h

Conservation of water surface elevation

When using the default adaptive reconstruction of variables, the Saint-Venant solver 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 η is conserved.

We start with the reconstruction of fine “wet” cells:

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

    double v = 0., η = 0.; // water surface elevation
    // 3x3 neighbourhood
    foreach_neighbor(1)
      if (h[] >= dry) {
	η += h[]*(zb[] + h[]);
	v += h[];
      }
    if (v > 0.)
      η /= v; // volume-averaged eta of neighbouring wet cells
    else

If none of the surrounding cells is wet, we assume a default sealevel at zero.

      η = 0.;

We then reconstruct the water depth in each child using η (of the parent cell i.e. a first-order interpolation in contrast to the wet case above) and zb of the child cells.

    // reconstruct water depth h from eta and zb
    foreach_child()
      h[] = max(0, η - zb[]);
  }
}

Cell restriction is simpler. We first compute the depth-weighted average of η over all the children…

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

… and use this in combination with zb (of the coarse cell) to compute the water depth h.

  if (v > 0.)
    h[] = max(0., η/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;
}
#else // Cartesian
void conserve_elevation (void) {}
#endif

“Radiation” boundary conditions

This can be used to implement open boundary conditions at low Froude numbers. The idea is to set the velocity normal to the boundary so that the water level relaxes towards its desired value (ref).

#define radiation(ref) (sqrt (G*max(h[],0.)) - sqrt(G*max((ref) - zb[], 0.)))

Tide gauges

An array of Gauge structures passed to output_gauges() will create a file (called name) for each gauge. Each time output_gauges() is called a line will be appended to the file. The line contains the time and the value of each scalar in list in the (wet) cell containing (x,y). The desc field can be filled with a longer description of the gauge.

typedef struct {
  char * name;
  double x, y;
  char * desc;
  FILE * fp;
} Gauge;

void output_gauges (Gauge * gauges, scalar * list)
{
  scalar * list1 = list_append (NULL, h);
  for (scalar s in list)
    list1 = list_append (list1, s);

  int n = 0;
  for (Gauge * g = gauges; g->name; g++, n++);
  coord a[n];
  n = 0;
  for (Gauge * g = gauges; g->name; g++, n++) {
    double xp = g->x, yp = g->y;
    unmap (&xp, &yp);
    a[n].x = xp, a[n].y = yp;
  }
  int len = list_len(list1);
  double v[n*len];
  interpolate_array (list1, a, n, v, false);

  if (pid() == 0) {
    n = 0;
    for (Gauge * g = gauges; g->name; g++) {
      if (!g->fp) {
	g->fp = fopen (g->name, "w");
	if (g->desc)
	  fprintf (g->fp, "%s\n", g->desc);
      }
      if (v[n] != nodata && v[n] > dry) {
	fprintf (g->fp, "%g", t);
	n++;
	for (scalar s in list)
	  fprintf (g->fp, " %g", v[n++]);
	fputc ('\n', g->fp);
	fflush (g->fp);
      }
      else
	n += len;
    }
  }

  free (list1);
}

Usage