#include "all-mach.h"
#include "tracer.h"
#include "elevation.h"
#include "gauges.h"

Include dependency graph for saint-venant-implicit.h:

Functions
void	event_defaults (void)
	Event: defaults (i = 0)

void	event_init (void)
	After user initialisation we set the initial pressure and apply boundary conditions.

void	event_stability (void)
	Event: stability (i++)

void	event_properties (void)
	The properties of the fluid are given by the "Saint-Venant equation of state" i.e.

void	eta_restriction (Point point, scalar s)
	On trees things are a bit more complicated due to the necessity to verify the "lake-at-rest" condition.

void	eta_prolongation (Point point, scalar s)
	Prolongation uses linear interpolation with a cell-centered gradient computed from simple finite-differences of wet cells.

void	event_acceleration (void)
	One can check that the prolongation values constructed using these functions verify \(\eta = constant\) for a lake-at-rest.

Variables
scalar	zb []
	In addition to the momentum \(\mathbf{q}=h\mathbf{u}\) defined by the all Mach solver, we will need the fluid depth h (i.e.

scalar	h []

scalar *	tracers = {q,h}
	Here we set the gradient functions for each tracer (as defined in the user-provided `tracers` list).

double	G = 1.

double	dry = 1e-4

scalar	rhoc2v []
	We need fields to store the (varying) state fields \(\rho c^2\), \(\alpha\) and the acceleration \(\mathbf{a}\).

vector	alphav []

vector	av []

double	CFLa = HUGE
	By default the timestep is only limited by the CFL condition on the advection velocity field.

double(*	gradient )(double, double, double) = minmod2
	The utility functions in elevation.h need to know which gradient we are using.

Function Documentation

◆ eta_prolongation()

void eta_prolongation	(	Point	point,
		scalar	s
	)

Prolongation uses linear interpolation with a cell-centered gradient computed from simple finite-differences of wet cells.

Definition at line 195 of file saint-venant-implicit.h.

References dimension, dry, g, h, s, vector::x, and x.

Referenced by event_acceleration().

Here is the caller graph for this function:

◆ eta_restriction()

void eta_restriction	(	Point	point,
		scalar	s
	)

On trees things are a bit more complicated due to the necessity to verify the "lake-at-rest" condition.

Topographic source term

The acceleration due to the topography is \(- g\nabla z_b\). On regular Cartesian grids we can simply do For a lake, the pressure gradient must balance the topographic source term i.e.

\[ \frac{1}{h}\nabla p = a = -g\nabla z_b \]

with \(a\) the acceleration. Using the equation of state and introducing

\[ \eta = z_b + h \]

the elevation of the free surface (constant for a lake at rest), we get

\[ \frac{1}{2h}\nabla gh^2 = -g\nabla (\eta - h) = g\nabla h \]

This identity is obviously verified mathematically, however it is not necessarily verified by the discrete gradient operator. In the case of Cartesian meshes it is simple to show that the naive discrete gradient operator we used above verifies this identity. For tree meshes this is not generally the case due to the prolongation operator used to fill ghost cells at refinement boundaries. Rather than trying to redefine the prolongation operator, we discretise the topographic source term as

\[ a = -g\nabla z_b = -g \nabla\eta + \frac{g}{2h}\nabla h^2 \]

This is strictly identical mathematically, however if we now use the same discrete operator to estimate \(\nabla p\) and \(\nabla h^2\), the discrete* lake-at-rest condition becomes

\[ \nabla\eta = 0 \]

which is much easier to verify.

To do so, we define the following restriction and prolongation operators for \(\eta\). The main idea is to avoid using the elevation of dry cells. Restriction is done by averaging the elevation of wet cells.

Definition at line 180 of file saint-venant-implicit.h.

References dry, h, n, nodata, s, sum, and x.

Referenced by event_acceleration().

Here is the caller graph for this function:

◆ event_acceleration()

void event_acceleration ( void )

One can check that the prolongation values constructed using these functions verify \(\eta = constant\) for a lake-at-rest.

Event: acceleration (i++)

To compute the acceleration due to topography, we first compute \(\eta\) and \(h^2\).

We then make sure that \(\eta\) uses our restriction/prolongation functions. \(h^2\) uses the same prolongation functions as \(p\) by default.

We then compute the acceleration as

\[ a = -g \nabla\eta + g \frac{\alpha}{2}\nabla h^2 \]

Definition at line 225 of file saint-venant-implicit.h.

References _i, alpha, av, eta, eta_prolongation(), eta_restriction(), G, h, sq(), vector::x, x, and zb.

Here is the call graph for this function:

◆ event_defaults()

void event_defaults ( void )

Event: defaults (i = 0)

Boundary conditions for VOF-advected tracers usually depend on boundary conditions for the VOF field.

Event: defaults (i = 0)

Initialisation

We set the default values.

Electrohydrodynamic stresses

The EHD force density, \(\mathbf{f}_e\), can be computed as the divergence of the Maxwell stress tensor \(\mathbf{M}\),

\[ M_{ij} = \varepsilon (E_i E_j - \frac{E^2}{2}\delta_{ij}) \]

where \(E_i\) is the \(i\)-component of the electric field, \(\mathbf{E}=-\nabla \phi\) and \(\delta_{ij}\) is the Kronecker delta.

We need to add the corresponding acceleration to the Navier–Stokes solver.

If the acceleration vector a (defined by the Navier–Stokes solver) is constant, we make it variable.

Event: defaults (i = 0)

Event: defaults (i = 0 )

Defaults

On trees we need to ensure conservation of the tracer when refining/coarsening.

Event: defaults (i = 0)

it is an acceleration. If necessary, we allocate a new vector field to store it.

Event: defaults (i = 0)

\[ \begin{aligned} 0 & = - \sum_k \nabla \cdot [\theta_H (hu)_k^{n + 1} + (1 - \theta_H) (hu)^n_k] \\ \frac{(hu)^{n + 1}_k - (hu)_k^n}{\Delta t} & = - \Delta tgh^{n + 1 / 2}_k (\theta_H \nabla \eta_r^{n + 1} + (1 - \theta_H) \nabla \eta_r^n) \end{aligned} \]

where \(\eta_r\) is the equivalent free-surface height (i.e. pressure) applied on the rigid lid.

Event: defaults (i = 0)

Setup

The \(w_k\) and \(\phi_k\) scalar fields are allocated and the \(w_k\) are added to the list of advected tracers.

Event: defaults (i = 0)

Initial conditions

Event: defaults (i = 0)

The default convention in Basilisk is no-flow through the boundaries of the domain, i.e. they are a streamline i.e. \(\psi=\)constant on the boundary. We set the default value for the CFL (the default in utils.h is 0.5). This is done once at the beginning of the simulation.

Event: defaults (i = 0)

We will need to add the acceleration term \(w^2/y\) in the evolution equation for \(u_y\). If the acceleration field is not allocated yet, we do so.

Event: defaults (i = 0)

We set limiting for q, h.

As well as default values.

On trees, we ensure that limiting is also applied to prolongation.

We setup the default display.

Definition at line 39 of file saint-venant-implicit.h.

References _i, a, alpha, alphav, av, display(), h, minmod2(), refine_linear(), rho, rhoc2, rhoc2v, s, and x.

Here is the call graph for this function:

◆ event_init()

void event_init ( void )

After user initialisation we set the initial pressure and apply boundary conditions.

The boundary conditions for \(\rho=h\) and \(\mathbf{q}\) are already applied by the all Mach solver.

Event: init (i = 0)

Definition at line 78 of file saint-venant-implicit.h.

References _i, G, h, p, sq(), and x.

Here is the call graph for this function:

◆ event_properties()

void event_properties ( void )

The properties of the fluid are given by the "Saint-Venant equation of state" i.e.

\(p = gh^2/2\), \(\rho c^2 = gh^2\).

Event: properties (i++)

The specific volume \(\alpha=1/\rho\) is constructed by averaging, taking into account dry states.

Definition at line 108 of file saint-venant-implicit.h.

References _i, alphav, dry, G, h, max, ps, rhoc2v, sq(), vector::x, x, and zb.

Here is the call graph for this function:

◆ event_stability()

void event_stability ( void )

Event: stability (i++)

We need to overload the stability event so that the CFL is taken into account (because we set stokes to true).

We first compute the minimum and maximum values of \(\alpha/f_m = 1/\rho\), as well as \(\Delta_{min}\).

The maximum timestep is set using the sum of surface tension coefficients.

Definition at line 93 of file saint-venant-implicit.h.

References CFLa, dry, dt, dtmax, G, h, HUGE, min, and x.

Variable Documentation

◆ alphav

vector alphav[]

Definition at line 36 of file saint-venant-implicit.h.

Referenced by event_defaults(), and event_properties().

◆ av

vector av[]

Definition at line 36 of file saint-venant-implicit.h.

Referenced by event_acceleration(), and event_defaults().

◆ CFLa

double CFLa = HUGE

By default the timestep is only limited by the CFL condition on the advection velocity field.

We add the option to define an "acoustic CFL number" which takes into account the speed of gravity waves.

Definition at line 90 of file saint-venant-implicit.h.

Referenced by event_stability().

◆ dry

double dry = 1e-4

Definition at line 29 of file saint-venant-implicit.h.

Referenced by eta_prolongation(), eta_restriction(), event_properties(), and event_stability().

◆ G

double G = 1.

Definition at line 28 of file saint-venant-implicit.h.

Referenced by event_acceleration(), event_init(), event_properties(), and event_stability().

◆ gradient

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

The utility functions in elevation.h need to know which gradient we are using.

Definition at line 263 of file saint-venant-implicit.h.

◆ h

scalar h[]

Definition at line 27 of file saint-venant-implicit.h.

Referenced by eta_prolongation(), eta_restriction(), event_acceleration(), event_defaults(), event_init(), event_properties(), and event_stability().

◆ rhoc2v

scalar rhoc2v[]

We need fields to store the (varying) state fields \(\rho c^2\), \(\alpha\) and the acceleration \(\mathbf{a}\).

Definition at line 35 of file saint-venant-implicit.h.

Referenced by event_defaults(), and event_properties().

◆ tracers

scalar * tracers = {q,h}

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 27 of file saint-venant-implicit.h.

◆ zb

scalar zb[]

In addition to the momentum \(\mathbf{q}=h\mathbf{u}\) defined by the all Mach solver, we will need the fluid depth h (i.e.

A semi-implicit Saint-Venant solver

Note that the multilayer solver provides the same functionality and should be prefered for most applications.

We solve the Saint-Venant equations semi-implicitly to lift the timestep restriction due to gravity waves. This is just a particular application of the "all Mach" semi-implicit solver. We will use the Bell-Collela-Glaz advection scheme to transport mass and momentum. the density \(\rho\)) and the topography zb. Both the momentum \(\mathbf{q}\) and mass \(h\) are advected tracers. The acceleration of gravity is G*. dry is the fluid depth below which a cell is considered "dry".

Definition at line 27 of file saint-venant-implicit.h.

Referenced by event_acceleration(), and event_properties().

Functions

Variables

Function Documentation

◆ eta_prolongation()

◆ eta_restriction()

Topographic source term

◆ event_acceleration()

◆ event_defaults()

Initialisation

Electrohydrodynamic stresses

Defaults

Setup

Initial conditions

◆ event_init()

◆ event_properties()

◆ event_stability()

Variable Documentation

◆ alphav

◆ av

◆ CFLa

◆ dry

◆ G

◆ gradient

◆ h

◆ rhoc2v

◆ tracers

Integral formulation for surface tension

◆ zb

A semi-implicit Saint-Venant solver