conservation.h

Go to the documentation of this file.

/** @file conservation.h
 */
/**
# Conservation of the volume of each layer
 
This module enforces conservation of the volume of each layer using a
spatially-constant, but varying in time, "upwelling" velocity within
each layer which compensates any volume loss. This can be interpreted
as a vertical volume/mass transfer between layers. 
 
This upwelling velocity is only applied when the thickness of a layer
is larger than `hmin/10`. The value of `hmin` is typically that set in
[entrainment.h](). */
 
extern double * hmin;
 
/**
The domain considered needs to be simply-connected in order to avoid
the transfer of mass between basins which are not connected. To
enforce this we use the functions in [/src/tag.h](). */
 
#include "tag.h"
 
/**
This stores the initial volume of each layer. */
 
static struct {
  double * sum;
} Conservation = {0};
 
/** @brief Event: init (i = 0) */
void event_init (void) 
{
 
  /**
  At initialisation we make sure that the domain is simply connected
  by retaining only the largest "oceanic" basin. We first define an
  indicator function which is zero on dry land and one in wet
  areas. */
  
  scalar d[];
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    double H = 0.;
    foreach_layer()
      H += h[];
    d[] = H > dry;
  }
 
  /**
  We remove cells which are only connected diagonally with their
  neighbors. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach */
    if ((d[] && !d[0,-1]) &&
  ((d[-1,-1] && !d[-1]) || (d[1,-1] && !d[1])))
      d[] = 0.;
 
  /**
  We only keep the largest basin. */
 
  remove_droplets (d, -1);
 
  /**
  We "dry out" all the other basins. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach */
    if (!d[]) {
      foreach_layer()
  h[] = 0.;
      zb[] = HUGE;
    }
 
  /**
  And finally we store the initial volume of each layer. */
 
  Conservation.sum = malloc (nl*sizeof (double));
  foreach_layer()
    Conservation.sum[_layer] = statsf (h).sum;
}
 
/**
At each timestep... */
 
/** @brief Event: viscous_term (i++) */
void event_viscous_term (void) 
{
  foreach_layer() {
 
    /**
    ... we compute the current volume of each layer and the area of
    cells which are "thick enough". */
    
    double area = 0., volume = 0.;
    foreach (reduction(+:area) reduction(+:volume)) {
      if (h[] > hmin[_layer]/10.)
  area += dv();
      volume += h[]*dv();
    }
 
    /**
    The ratio of these two quantities gives the vertical displacement,
    which is then applied in each cell which is "thick enough". */
 
    assert (area > 0.);
    double dh = (Conservation.sum[_layer] - volume)/area;
    assert (dh > - hmin[_layer]/10.);
    for (int _i = 0; _i < _N; _i++) /* foreach */
      if (h[] > hmin[_layer]/10.)
  h[] += dh;
  }
}
 
/**
We free memory at the end of the simulation to avoid memory leaks. */
 
/** @brief Event: cleanup (t = end) */
void event_cleanup (void) 
{
  free (Conservation.sum);
}