hydro.h File Reference

#include "utils.h"
#include "run.h"
#include "diffusion.h"
#include "elevation.h"
#include "gauges.h"

Include dependency graph for hydro.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Data Structures
struct	Flux
	A NULL-terminated array of Flux structures passed to output_fluxes()* will create a file (called name) for each flux. More...

Macros
#define	BGHOSTS   2
#define	LAYERS   1
#define	gmetric(i)   (2.*fm.x[i]/(cm[i] + cm[i-1]))
	The macro below can be overloaded to define the barotropic acceleration.
#define	a_baro(eta, i)   (G*gmetric(i)*(eta[i-1] - eta[i])/Delta)
#define	slope_limited(dz)
#define	hpg(pg, phi, i, code)
#define	radiation(ref)   _radiation(point, ref, _s)
#define	conserve_elevation()   conserve_layered_elevation()

Functions
static void	refine_eta (Point point, scalar eta)
static void	restriction_eta (Point point, scalar eta)
void	event_defaults0 (void)
	The allocation of fields for each layer is split in two parts, because we need to make sure that the layer thicknesses and \(\eta\) are allocated first (in the case where other layer fields are added, for example for the non-hydrostatic extension).
void	event_defaults (void)
	Other fields, such as \(\mathbf{u}_k\) here, are added by this event.
void	event_init (void)
	After user initialisation, we define the free-surface height \(\eta\).
event	set_dtmax (i++, last) dtmax
void	event_face_fields (void)
	Event: face_fields (i++, last)
trace void	advect (scalar *tracers, vector hu, vector hf, double dt)
event	half_advection (i++, last)
	This is where the 'two-step advection' of the implicit scheme plugs itself (nothing is done for the explicit scheme).
event	acceleration (i++, last)
	Vertical diffusion (including viscosity) is added by this code.
void	event_pressure (void)
	Event: pressure (i++, last)
void	event_update_eta (void)
	Finally the free-surface height \(\eta\) is updated.
event	remap (i++, last)
event	adapt (i++, last)
void	event_cleanup (void)
	Event: cleanup (t = end, last)
void	vertical_velocity (scalar w, vector hu, vector hf)
double	_radiation (Point point, double ref, scalar s)
static void	refine_layered_elevation (Point point, scalar h)
	But we need to re-define the water depth refinement function.
void	conserve_layered_elevation (void)
	We overload the conserve_elevation() function.
double	segment_flux (coord segment[2], double *flux, scalar h, vector u)
void	output_fluxes (Flux *fluxes, scalar h, vector u)

Variables
scalar	zb []
scalar	eta
scalar	h
vector	u
double	G = 1.
double	CFL_H = 1e40
double	dry = 1e-12
double(*	gradient )(double, double, double) = minmod2
scalar *	tracers = NULL
	Here we set the gradient functions for each tracer (as defined in the user-provided tracers list).
scalar	eta_r
bool	linearised = false
double	dtmax
	The maximum timestep dtmax can be used to impose additional stability conditions.
static bool	hydrostatic = true
vector	hu
vector	hf
vector	ha
double	max_slope = 0.577350269189626 [0]

Macro Definition Documentation

a_baro

#define a_baro	(		eta,
			i
	)		(G*gmetric(i)*(eta[i-1] - eta[i])/Delta)

Definition at line 197 of file hydro.h.

BGHOSTS

#define BGHOSTS   2

The hydrostatic multilayer solver for free-surface flows

The theoretical basis and main algorithms for this solver are described in Popinet, 2020. Note however that this version of the solver is more recent and may not match the details of Popinet, 2020.

The system of \(n\) layers of incompressible fluid pictured below is modelled as the set of (hydrostatic) equations:

\[ \begin{aligned} \partial_t h_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \right)_k & = 0,\\ \partial_t \left( h \mathbf{u} \right)_k + \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \mathbf{u} \right)_k & = - gh_k \mathbf{{\nabla}} (\eta) \end{aligned} \]

with \(\mathbf{u}_k\) the velocity vector, \(h_k\) the layer thickness, \(g\) the acceleration of gravity and \(\eta\) the free-surface height. The non-hydrostatic pressure \(\phi_k\) and vertical velocity \(w_k\) can be added using the non-hydrostatic extension. See also the [general introduction](README).

Fields and parameters

The zb and eta fields define the bathymetry and free-surface height. The h and u fields define the layer thicknesses and velocities.

The acceleration of gravity is G and dry controls the minimum layer thickness. The hydrostatic CFL criterion is defined by CFL_H. It is set to a very large value by default, but will be tuned either by the user or by the default solver settings (typically depending on whether time integration is explicit or implicit).

The gradient pointer gives the function used for limiting.

tracers is a list of tracers for each layer. By default it contains only the components of velocity (unless linearised is set to true in which case the tracers list is empty).

Definition at line 46 of file hydro.h.

conserve_elevation

#define conserve_elevation

(

void

)

conserve_layered_elevation()

Definition at line 658 of file hydro.h.

gmetric

#define gmetric

(

i

)

(2.*fm.x[i]/(cm[i] + cm[i-1]))

The macro below can be overloaded to define the barotropic acceleration.

By default it is just the slope of the free-surface times gravity.

Definition at line 195 of file hydro.h.

hpg

#define hpg	(		pg,
			phi,
			i,
			code
	)

Value:

  do {            \
  double dz = zb[i] - zb[i-1];            \
  foreach_layer() {             \
    double pg = 0.;             \
    if (h[i] + h[i-1] > dry) {            \
      double s = Delta*slope_limited(dz/Delta);       \
      pg = (h[i-1] - s)*phi[i-1] - (h[i] + s)*phi[i];                   \
      if (point.l < nl - 1) {           \
  double s = Delta*slope_limited((dz + h[i] - h[i-1])/Delta);     \
  pg += (h[i-1] + s)*phi[i-1,0,1] - (h[i] - s)*phi[i,0,1];        \
      }                 \
      pg *= gmetric(i)*hf.x[i]/(Delta*(h[i] + h[i-1]));     \
    } code;               \
    dz += h[i] - h[i-1];            \
  }                 \
} while (0)

Definition at line 525 of file hydro.h.

LAYERS

#define LAYERS   1

Definition at line 47 of file hydro.h.

radiation

#define radiation

(

ref

)

_radiation(point, ref, _s)

Definition at line 593 of file hydro.h.

slope_limited

#define slope_limited

(

dz

)

Value:

         (fabs(dz) < max_slope ? (dz) : \
         ((dz) > 0. ? max_slope : - max_slope))

Definition at line 522 of file hydro.h.

Function Documentation

_radiation()

double _radiation	(	Point	point,
		double	ref,
		scalar	s
	)

"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).

Definition at line 581 of file hydro.h.

References dry, foreach_layer, G, h, max, x, and zb.

acceleration()

event acceleration	(	i++	,
		last
	)

Vertical diffusion (including viscosity) is added by this code.

Accelerations

Acceleration terms are added here. In the simplest case, this is only the pressure gradient due to the free-surface slope, as computed in face_fields.

adapt()

event adapt	(	i++	,
		last
	)

advect()

trace void advect	(	scalar *	tracers,
		vector	hu,
		vector	hf,
		double	dt
	)

Advection and diffusion

The function below approximates the advection terms using estimates of the face fluxes \(h\mathbf{u}\) and face heights \(h_f\).

The fluxes are first limited according to the CFL condition to ensure strict positivity of the layer heights. This step is necessary due to the approximate estimation of the CFL condition in the timestep calculation above.

In the case where the flux is limited, and for multiple layers, an attempt is made to conserve the total barotropic flux by merging the flux difference with the flux in the layer just above. This allows to maintain an accurate evolution for the free-surface height \(\eta\).

We compute the flux _{i+1/2,k}$ for each tracer \(s\), using a variant of the BCG scheme.

We compute ^\star_i = (hs)^n_i + \Delta t [(shu)_{i+1/2} -(shu)_{i-1/2}]/\Delta$.

We then obtain \(h^{n+1}\) and \(s^{n+1}\) using

\[ \begin{aligned} h_i^{n+1} & = h_i^n + \Delta t \frac{(hu)_{i+1/2} - (hu)_{i-1/2}}{\Delta},\\ s_i^{n+1} & = \frac{(hs)^\star_i}{h_i^{n+1}} \end{aligned} \]

Definition at line 306 of file hydro.h.

References _i, _layer, a, CFL, cm, dimension, dry, dt, f, flux, foreach_layer, g, h, hf, hu, i, LINENO, max, nl, point, s, sign(), t, un, vector::x, x, vector::y, and y.

Referenced by event_half_advection(), and event_pressure().

Here is the call graph for this function:

Here is the caller graph for this function:

conserve_layered_elevation()

void conserve_layered_elevation

(

void

)

We overload the conserve_elevation() function.

Definition at line 651 of file hydro.h.

References h, prolongation_elevation(), refine_layered_elevation(), restriction_elevation(), set_prolongation(), and set_restriction().

Here is the call graph for this function:

event_cleanup()

void event_cleanup

(

void

)

Event: cleanup (t = end, last)

Cleanup

The fields and lists allocated in `defaults()` above must be freed at the end of the run.

Definition at line 503 of file hydro.h.

References eta, eta_r, free(), h, tracers, u, and x.

Here is the call graph for this function:

event_defaults()

void event_defaults

(

void

)

Other fields, such as \(\mathbf{u}_k\) here, are added by this event.

Event: defaults (i = 0)

The (velocity) CFL is limited by the unsplit advection scheme, so is dependent on the dimension. The (gravity wave) CFL is set to 1/2 (if not already set by the user).

The gradient and prolongation/restriction functions are set for all tracer fields.

We setup the default display.

Definition at line 131 of file hydro.h.

References CFL, CFL_H, dimension, display(), gradient, linearised, list_append(), nl, refine_linear(), reset, restriction_volume_average(), s, set_prolongation(), set_restriction(), tracers, u, vector::x, and x.

Here is the call graph for this function:

event_defaults0()

void event_defaults0

(

void

)

The allocation of fields for each layer is split in two parts, because we need to make sure that the layer thicknesses and \(\eta\) are allocated first (in the case where other layer fields are added, for example for the non-hydrostatic extension).

Event: defaults0 (i = 0)

We set the proper gradient and refinement/restriction functions.

Definition at line 95 of file hydro.h.

References assert, eta, eta_r, gradient, h, list_copy(), nl, refine_eta(), refine_linear(), reset, restriction_eta(), restriction_volume_average(), set_prolongation(), set_restriction(), and zb.

Here is the call graph for this function:

event_face_fields()

void event_face_fields

(

void

)

Event: face_fields (i++, last)

The (CFL-limited) timestep is also computed by this function. A difficulty is that the prediction step below also requires an estimated timestep (the pdt variable below). The timestep at the previous iteration is used as estimate.

The face velocity is computed as the height-weighted average of the cell velocities.

If the left or central cell are dry, we consider a "step-like" bathymetry and define the face height as the water level above the step.

In the default case, the face height is computed using a variant of the BCG scheme.

The maximum velocity is stored and the flux and height-weighted accelerations are computed.

The maximum timestep is computed using the total depth H and the advection and gravity wave CFL criteria. The gravity wave speed takes dispersion into account in the non-hydrostatic case.

The timestep is computed, taking into account the timing of events, and also stored in pdt (see comment above).

Definition at line 212 of file hydro.h.

References _i, a, a_baro, c, CFL, CFL_H, cm, dry, dt, dtmax, dtnext(), eta_r, fm, foreach_layer, g, G, h, ha, hf, hu, hydrostatic, i, min, nl, sign(), u, un, vector::x, x, and zb.

Here is the call graph for this function:

event_init()

void event_init

(

void

)

After user initialisation, we define the free-surface height \(\eta\).

Event: init (i = 0)

Definition at line 173 of file hydro.h.

References _i, eta, foreach_layer, h, x, and zb.

event_pressure()

void event_pressure

(

void

)

Event: pressure (i++, last)

The acceleration is applied to the face fluxes...

... and to the centered velocity field, using height-weighting.

The resulting fluxes are used to advect both tracers and layer heights.

Definition at line 434 of file hydro.h.

References _i, advect(), cm, dimension, dry, dt, fm, foreach_layer, ha, hf, hu, tracers, u, vector::x, x, and vector::y.

Here is the call graph for this function:

event_update_eta()

void event_update_eta

(

void

)

Finally the free-surface height \(\eta\) is updated.

Event: update_eta (i++, last)

Definition at line 474 of file hydro.h.

References _i, eta, foreach_layer, h, hf, hu, x, and zb.

half_advection()

event half_advection	(	i++	,
		last
	)

This is where the 'two-step advection' of the implicit scheme plugs itself (nothing is done for the explicit scheme).

output_fluxes()

void output_fluxes	(	Flux *	fluxes,
		scalar	h,
		vector	u
	)

Definition at line 726 of file hydro.h.

References f, fflush(), flux, h, i, nl, segment_flux(), t, u, and x.

Here is the call graph for this function:

refine_eta()

static void refine_eta ( Point  point, scalar  eta  )

static

Setup

We first define refinement and restriction functions for the free-surface height eta which is just the sum of all layer thicknesses and of bathymetry.

Definition at line 71 of file hydro.h.

References eta, foreach_layer, h, x, and zb.

Referenced by event_defaults0().

Here is the caller graph for this function:

refine_layered_elevation()

static void refine_layered_elevation ( Point  point, scalar  h  )

static

But we need to re-define the water depth refinement function.

Conservation of water surface elevation

We re-use some generic functions.

We first check whether we are dealing with "shallow cells".

If we do, refined cells are just set to the default sea level.

Otherwise, we use the surface elevation of the parent cells to reconstruct the water depth of the children cells.

Definition at line 606 of file hydro.h.

References default_sea_level, dimension, eta, g, gradient, h, max, vector::x, x, and zb.

Referenced by conserve_layered_elevation().

Here is the caller graph for this function:

remap()

event remap	(	i++	,
		last
	)

Vertical remapping and mesh adaptation

Vertical remapping is applied here if necessary.

restriction_eta()

static void restriction_eta ( Point  point, scalar  eta  )

static

Definition at line 80 of file hydro.h.

References eta, foreach_layer, h, and zb.

Referenced by event_defaults0(), and event_pressure().

Here is the caller graph for this function:

segment_flux()

double segment_flux	(	coord	segment[2],
		double *	flux,
		scalar	h,
		vector	u
	)

Fluxes through sections

These functions are typically used to compute fluxes (i.e. flow rates) through cross-sections defined by two endpoints (i.e. segments). Note that the orientation of the segment is taken into account when computing the flux i.e the positive normal direction to the segment is to the "left" when looking from the start to the end.

This can be expressed mathematically as:

\[ \text{flux}[k] = \int_A^B h_k\mathbf{u}_k\cdot\mathbf{n}dl \]

with \(A\) and \(B\) the endpoints of the segment, \(k\) the layer index, \(\mathbf{n}\) the oriented segment unit normal and \(dl\) the elementary length. The function returns the sum (over \(k\)) of all the fluxes.

Definition at line 683 of file hydro.h.

References a, dimension, dp, flux, fm, foreach_layer, foreach_segment(), h, i, interpolate_linear(), l, m(), nl, normalize(), p, point, sq(), u, vector::x, x, vector::y, and coord::y.

Referenced by output_fluxes().

Here is the call graph for this function:

Here is the caller graph for this function:

set_dtmax()

event set_dtmax	(	i++	,
		last
	)

vertical_velocity()

void vertical_velocity	(	scalar	w,
		vector	hu,
		vector	hf
	)

Hydrostatic vertical velocity

For the hydrostatic solver, the vertical velocity is not defined by default since it is usually not required. The function below can be applied to compute it using the (diagnostic) incompressibility condition

\[ \mathbf{{\nabla}} \cdot \left( h \mathbf{u} \right)_k + \left[ w - \mathbf{u} \cdot \mathbf{{\nabla}} (z) \right]_k = 0 \]

Definition at line 554 of file hydro.h.

References _i, dimension, dry, foreach_layer, h, hf, hu, point, w, vector::x, x, and zb.

Variable Documentation

CFL_H

double CFL_H = 1e40

Definition at line 52 of file hydro.h.

Referenced by event_defaults(), event_defaults0(), and event_face_fields().

dry

double dry = 1e-12

Definition at line 56 of file hydro.h.

Referenced by _radiation(), advect(), event_face_fields(), event_pressure(), and vertical_velocity().

dtmax

double dtmax

The maximum timestep dtmax can be used to impose additional stability conditions.

Definition at line 187 of file hydro.h.

Referenced by event_face_fields().

eta

scalar eta

Definition at line 50 of file hydro.h.

Referenced by box_matrix(), event_acceleration(), event_cleanup(), event_defaults0(), event_face_fields(), event_init(), event_pressure(), event_tracer_advection(), event_update_eta(), event_viscous_term(), fault(), fenep(), geostrophic_velocity(), rectangular_source(), refine_elevation(), refine_eta(), refine_layered_elevation(), relax_hydro(), relax_nh(), residual_hydro(), residual_nh(), restriction_elevation(), restriction_eta(), turbulence_fk_craig(), and update_green_naghdi().

eta_r

scalar eta_r

Definition at line 60 of file hydro.h.

Referenced by event_acceleration(), event_cleanup(), event_defaults(), event_defaults0(), event_face_fields(), and event_pressure().

G

double G = 1.

Definition at line 52 of file hydro.h.

Referenced by _radiation(), and event_face_fields().

gradient

double(* gradient) (double, double, double)	(	double	,
		double	,
		double
	)		= minmod2

Definition at line 58 of file hydro.h.

Referenced by event_defaults(), event_defaults0(), and refine_layered_elevation().

h

scalar h

Definition at line 50 of file hydro.h.

Referenced by _radiation(), advect(), conserve_layered_elevation(), event_cleanup(), event_defaults0(), event_face_fields(), event_init(), event_update_eta(), output_fluxes(), refine_eta(), refine_layered_elevation(), restriction_eta(), segment_flux(), and vertical_velocity().

ha

vector ha

Definition at line 209 of file hydro.h.

Referenced by de(), event_acceleration(), event_face_fields(), event_pressure(), event_viscous_term(), and sum_add_sum().

hf

vector hf

Definition at line 209 of file hydro.h.

Referenced by advect(), event_acceleration(), event_face_fields(), event_half_advection(), event_pressure(), event_update_eta(), event_viscous_term(), horizontal_diffusion(), implicit_horizontal_diffusion(), and vertical_velocity().

hu

vector hu

Definition at line 209 of file hydro.h.

Referenced by advect(), event_acceleration(), event_face_fields(), event_half_advection(), event_pressure(), event_update_eta(), and vertical_velocity().

hydrostatic

bool hydrostatic = true

static

Computation of face values

At each timestep, temporary face fields are defined for the fluxes $(h \mathbf{u})^{n+1/2}$, face height \(h_f^{n+1/2}\) and height-weighted face accelerations ^{n+1/2}$.

Definition at line 208 of file hydro.h.

Referenced by event_defaults(), and event_face_fields().

linearised

bool linearised = false

Definition at line 61 of file hydro.h.

Referenced by event_defaults().

max_slope

double max_slope = 0.577350269189626 [0]

Horizontal pressure gradient

The macro below computes the horizontal pressure gradient

\[ pg_k = - \mathbf{{\nabla}} (h \phi)_k + \left[ \phi \mathbf{{\nabla}} z \right]_k \]

on the faces of each layer. The slope of the layer interfaces \(\mathbf{{\nabla}} z_{k+1/2}\) in the second-term is bounded by max_slope (by default 30 degrees).

Definition at line 521 of file hydro.h.

tracers

scalar* tracers = NULL

Here we set the gradient functions for each tracer (as defined in the user-provided tracers list).

Integral formulation for surface tension

See Al Saud et al., 2018 and Popinet & Zaleski, 1999 for details.

The surface tension field \(\sigma\) will be associated to each levelset tracer. This is done easily by adding the following field attributes.

Definition at line 60 of file hydro.h.

Referenced by event_cleanup(), event_defaults(), and event_pressure().

u

vector u

Definition at line 51 of file hydro.h.

Referenced by event_cleanup(), event_defaults(), event_face_fields(), event_pressure(), output_fluxes(), and segment_flux().

zb

scalar zb[]

Definition at line 50 of file hydro.h.

Referenced by _radiation(), event_defaults0(), event_face_fields(), event_init(), event_update_eta(), refine_eta(), refine_layered_elevation(), restriction_eta(), and vertical_velocity().