#include "run.h"
#include "timestep.h"
#include "bcg.h"
#include "viscosity.h"

Include dependency graph for centered.h:

Macros
#define	neumann_pressure(i) (a.n[i]*fm.n[i]/alpha.n[i])

Functions
void	event_defaults (void)
	For embedded boundaries on trees, we need to define the pressure gradient for prolongation of pressure close to embedded boundaries.

void	event_default_display (void) display("squares (color
	We had some objects to display by default.

void	event_init (void)
	Event: init (i = 0)

event	set_dtmax (i++, last) dtmax

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

event	vof (i++, last)
	If we are using VOF or diffuse tracers, we need to advance them (to time \(t+\Delta t/2\)) here.

event	tracer_advection (i++, last)

event	tracer_diffusion (i++, last)

event	properties (i++, last)
	The fluid properties such as specific volume (fields \(\alpha\) and \(\alpha_c\)) or dynamic viscosity (face field \(\mu_f\)) – at time \(t+\Delta t/2\) – can be defined by overloading this event.

void	prediction ()

void	event_advection_term (void)
	Event: advection_term (i++,last)

static void	correction (double dt)

void	event_viscous_term (void)
	The viscous term is computed implicitly.

void	event_acceleration (void)
	Event: acceleration (i++,last)

void	centered_gradient (scalar p, vector g)

void	event_projection (void)
	To get the pressure field at time \(t + \Delta t\) we project the face velocity field (which will also be used for tracer advection at the next timestep).

event	end_timestep (i++, last)
	Some derived solvers need to hook themselves at the end of the timestep.

void	event_adapt (void)
	Event: adapt (i++,last)

Variables
scalar	p [] = neumann ((a.n[ ghost ]*fm.n[ ghost ]/alpha.n[ ghost ]))
	The primary variables are the centered pressure field \(p\) and the centered velocity field \(\mathbf{u}\).

vector	uf []
	We allocate the (face) velocity field.

vector	g []

scalar	pf []

const vector	mu = zerof
	In the case of variable density, the user will need to define both the face and centered specific volume fields ( \(\alpha\) and \(\alpha_c\) respectively) i.e.

const vector	a = zerof

const vector	alpha = unityf

const scalar	rho = unity

mgstats	mgp = {0}

mgstats	mgpf = {0}

mgstats	mgu = {0}

bool	stokes = false

void	spread = -1)

double	dtmax
	After user initialisation, we initialise the face velocity and fluid properties.

Macro Definition Documentation

◆ neumann_pressure

#define neumann_pressure ( i ) (a.n[i]*fm.n[i]/alpha.n[i])

Boundary conditions

For the default symmetric boundary conditions, we need to ensure that the normal component of the velocity is zero after projection. This means that, at the boundary, the acceleration \(\mathbf{a}\) must be balanced by the pressure gradient. Taking care of boundary orientation and staggering of \(\mathbf{a}\), this can be written

Definition at line 90 of file centered.h.

Function Documentation

◆ centered_gradient()

void centered_gradient	(	scalar	p,
		vector	g
	)

Approximate projection

This function constructs the centered pressure gradient and acceleration field g using the face-centered acceleration field a and the cell-centered pressure field p.

We first compute a face field \(\mathbf{g}_f\) combining both acceleration and pressure gradient.

We average these face values to obtain the centered, combined acceleration and pressure gradient field.

Definition at line 411 of file centered.h.

References _i, a, alpha, dimension, fm, g, p, SEPS, vector::x, and x.

Referenced by event_end_timestep(), and event_projection().

Here is the caller graph for this function:

◆ correction()

static void correction ( double dt )

static

Viscous term

We first define a function which adds the pressure gradient and acceleration terms.

Definition at line 345 of file centered.h.

References _i, dimension, dt, g, u, vector::x, and x.

Referenced by event_acceleration(), event_projection(), and event_viscous_term().

Here is the caller graph for this function:

◆ end_timestep()

event end_timestep	(	i++	,
		last
	)

Some derived solvers need to hook themselves at the end of the timestep.

◆ event_acceleration()

void event_acceleration ( void )

Event: acceleration (i++,last)

Acceleration term

The acceleration term \(\mathbf{a}\) needs careful treatment as many equilibrium solutions depend on exact balance between the acceleration term and the pressure gradient: for example Laplace's balance for surface tension or hydrostatic pressure in the presence of gravity.

To ensure a consistent discretisation, the acceleration term is defined on faces as are pressure gradients and the centered combined acceleration and pressure gradient term \(\mathbf{g}\) is obtained by averaging.

The (provisionary) face velocity field at time \(t+\Delta t\) is obtained by interpolation from the centered velocity field. The acceleration term is added.

Definition at line 397 of file centered.h.

References _i, a, dt, face_value, fm, u, uf, vector::x, and x.

◆ event_adapt()

void event_adapt ( void )

Event: adapt (i++,last)

Adaptivity

After mesh adaptation fluid properties need to be updated. When using embedded boundaries the fluid fractions and face fluxes need to be checked for inconsistencies.

Definition at line 464 of file centered.h.

References _i, cs, event(), fractions_cleanup(), fs, uf, vector::x, and x.

Here is the call graph for this function:

◆ event_advection_term()

void event_advection_term ( void )

Event: advection_term (i++,last)

Advection term

We predict the face velocity field \(\mathbf{u}_f\) at time \(t+\Delta t/2\) then project it to make it divergence-free. We can then use it to compute the velocity advection term, using the standard Bell-Collela-Glaz advection scheme for each component of the velocity field.

Definition at line 330 of file centered.h.

References advection(), alpha, dt, g, mgpf, mgstats::nrelax, pf, prediction(), project(), stokes, u, and uf.

Here is the call graph for this function:

◆ event_default_display()

void event_default_display ( void )

We had some objects to display by default.

Event: default_display (i = 0)

◆ event_defaults()

void event_defaults ( void )

For embedded boundaries on trees, we need to define the pressure gradient for prolongation of pressure close to embedded boundaries.

Event: defaults (i = 0)

Initial conditions

Event: defaults (i = 0)

We reset the multigrid parameters to their default values.

The pressures are never dumped.

The default density field is set to unity (times the metric and the solid factors).

On trees, refinement of the face-centered velocity field needs to preserve the divergence-free condition.

When using embedded boundaries, the restriction and prolongation operators need to take the boundary into account.

We set the dimensions of the velocity field.

Definition at line 129 of file centered.h.

References _i, alpha, alphav, CFL, cm, cs, dimension, fm, scalar::i, list_add(), mgp, mgpf, mgu, p, pf, refine_embed_linear(), refine_face(), refine_face_solenoidal(), restriction_embed_linear(), rho, s, uf, unityf, vector::x, and x.

Here is the call graph for this function:

◆ event_init()

void event_init ( void )

Event: init (i = 0)

We update fluid properties.

We set the initial timestep (this is useful only when restoring from a previous run).

Definition at line 208 of file centered.h.

References _i, DT, dtmax, event(), face_value, fm, u, uf, vector::x, and x.

Here is the call graph for this function:

◆ event_projection()

void event_projection ( void )

To get the pressure field at time \(t + \Delta t\) we project the face velocity field (which will also be used for tracer advection at the next timestep).

Then compute the centered gradient field g.

Event: projection (i++,last)

We add the gradient field g to the centered velocity field.

Definition at line 438 of file centered.h.

References alpha, centered_gradient(), correction(), dt, g, mgp, mgstats::nrelax, p, project(), and uf.

Here is the call graph for this function:

◆ event_stability()

void event_stability ( void )

Event: stability (i++,last)

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 237 of file centered.h.

References dt, dtmax, dtnext(), stokes, timestep(), and uf.

Here is the call graph for this function:

◆ event_viscous_term()

void event_viscous_term ( void )

The viscous term is computed implicitly.

We first add the pressure gradient and acceleration terms, as computed at time \(t\), then call the implicit viscosity solver. We then remove the acceleration and pressure gradient terms as they will be replaced by their values at time \(t+\Delta t\).

Event: viscous_term (i++,last)

We reset the acceleration field (if it is not a constant).

Definition at line 360 of file centered.h.

References _i, a, correction(), dt, mgu, mu, mgstats::nrelax, rho, u, viscosity(), vector::x, and x.

Here is the call graph for this function:

◆ prediction()

void prediction ( )

Predicted face velocity field

For second-order in time integration of the velocity advection term \(\nabla\cdot(\mathbf{u}\otimes\mathbf{u})\), we need to define the face velocity field \(\mathbf{u}_f\) at time \(t+\Delta t/2\). We use a version of the Bell-Collela-Glaz advection scheme and the pressure gradient and acceleration terms at time \(t\) (stored in vector \(\mathbf{g}\)).

Definition at line 268 of file centered.h.

References _i, dimension, dt, fm, fs, g, i, s, sign(), u, uf, un, vector::x, x, and vector::y.

Referenced by event_advection_term().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ properties()

event properties	(	i++	,
		last
	)

The fluid properties such as specific volume (fields \(\alpha\) and \(\alpha_c\)) or dynamic viscosity (face field \(\mu_f\)) – at time \(t+\Delta t/2\) – can be defined by overloading this event.

◆ set_dtmax()

event set_dtmax	(	i++	,
		last
	)

Time integration

The timestep for this iteration is controlled by the CFL condition, applied to the face centered velocity field \(\mathbf{u}_f\); and the timing of upcoming events.

◆ tracer_advection()

event tracer_advection	(	i++	,
		last
	)

◆ tracer_diffusion()

event tracer_diffusion	(	i++	,
		last
	)

◆ vof()

event vof	(	i++	,
		last
	)

If we are using VOF or diffuse tracers, we need to advance them (to time \(t+\Delta t/2\)) here.

Note that this assumes that tracer fields are defined at time \(t-\Delta t/2\) i.e. are lagging the velocity/pressure fields by half a timestep.

Variable Documentation

◆ a

const vector a = zerof

Definition at line 72 of file centered.h.

Referenced by centered_gradient(), event_acceleration(), and event_viscous_term().

◆ alpha

const vector alpha = unityf

Definition at line 72 of file centered.h.

Referenced by centered_gradient(), event_advection_term(), event_defaults(), and event_projection().

◆ dtmax

double dtmax

After user initialisation, we initialise the face velocity and fluid properties.

Definition at line 205 of file centered.h.

Referenced by event_init(), and event_stability().

◆ g

vector g[]

Definition at line 47 of file centered.h.

Referenced by centered_gradient(), correction(), event_advection_term(), event_projection(), and prediction().

◆ mgp

mgstats mgp = {0}

Definition at line 74 of file centered.h.

Referenced by event_defaults(), and event_projection().

◆ mgpf

mgstats mgpf = {0}

Definition at line 74 of file centered.h.

Referenced by event_advection_term(), event_defaults(), and event_perfs().

◆ mgu

mgstats mgu = {0}

Definition at line 74 of file centered.h.

Referenced by event_defaults(), and event_viscous_term().

◆ mu

const vector mu = zerof

In the case of variable density, the user will need to define both the face and centered specific volume fields ( \(\alpha\) and \(\alpha_c\) respectively) i.e.

\(1/\rho\). If not specified by the user, these fields are set to one i.e. the density is unity.

Viscosity is set by defining the face dynamic viscosity \(\mu\); default is zero.

The face field \(\mathbf{a}\) defines the acceleration term; default is zero.

The statistics for the (multigrid) solution of the pressure Poisson problems and implicit viscosity are stored in mgp, mgpf, mgu respectively.

If stokes is set to true, the velocity advection term \(\nabla\cdot(\mathbf{u}\otimes\mathbf{u})\) is omitted. This is a reference to Stokes flows for which inertia is negligible compared to viscosity.

Definition at line 72 of file centered.h.

Referenced by event_viscous_term().

◆ p

p[bottom] = neumann ((a.n[ ghost ]*fm.n[ ghost ]/alpha.n[ ghost ]))

The primary variables are the centered pressure field \(p\) and the centered velocity field \(\mathbf{u}\).

Incompressible Navier–Stokes solver (centered formulation)

We wish to approximate numerically the incompressible, variable-density Navier–Stokes equations

\[ \partial_t\mathbf{u}+\nabla\cdot(\mathbf{u}\otimes\mathbf{u}) = \frac{1}{\rho}\left[-\nabla p + \nabla\cdot(2\mu\mathbf{D})\right] + \mathbf{a} \]

\[ \nabla\cdot\mathbf{u} = 0 \]

with the deformation tensor \(\mathbf{D}=[\nabla\mathbf{u} + (\nabla\mathbf{u})^T]/2\).

The scheme implemented here is close to that used in Gerris (Popinet, 2003, Popinet, 2009, Lagrée et al, 2011).

We will use the generic time loop, a CFL-limited timestep, the Bell-Collela-Glaz advection scheme and the implicit viscosity solver. If embedded boundaries are used, a different scheme is used for viscosity. The centered vector field \(\mathbf{g}\) will contain pressure gradients and acceleration terms.

We will also need an auxilliary face velocity field \(\mathbf{u}_f\) and the associated centered pressure field \(p_f\).

Definition at line 46 of file centered.h.

Referenced by centered_gradient(), event_defaults(), and event_projection().

◆ pf

scalar pf[]

Definition at line 48 of file centered.h.

Referenced by event_advection_term(), and event_defaults().

◆ rho

const scalar rho = unity

Definition at line 73 of file centered.h.

Referenced by event_defaults(), and event_viscous_term().

◆ spread

void spread = -1)

Definition at line 199 of file centered.h.

Referenced by if(), if(), if(), and if().

◆ stokes

bool stokes = false

Definition at line 75 of file centered.h.

Referenced by event_advection_term(), event_defaults(), and event_stability().

◆ uf

vector uf[]

We allocate the (face) velocity field.

An advection solver

We wish to solve the advection equations

\[ \partial_tf_i+\mathbf{u}\cdot\nabla f_i=0 \]

where \(\mathbf{u}\) is the velocity field and \(f_i\) are a list of passive tracers. This can be done with a flux-based advection scheme such as the 2nd-order, unsplit, upwind scheme of Bell-Collela-Glaz, 1989.

The main time loop is defined in [run.h](). A stable timestep needs to respect the CFL condition. For compatibility with the other solvers, we allocate it as uf and define an alias. The gradient function is used to set the type of slope-limiting required. The default is to not use any limiting (i.e. a purely centered slope estimation).

Definition at line 47 of file centered.h.

Referenced by event_acceleration(), event_adapt(), event_advection_term(), event_defaults(), event_init(), event_projection(), event_stability(), event_tracer_advection(), for(), and prediction().

Macros

Functions

Variables

Macro Definition Documentation

◆ neumann_pressure

Boundary conditions

Function Documentation

◆ centered_gradient()

Approximate projection

◆ correction()

Viscous term

◆ end_timestep()

◆ event_acceleration()

Acceleration term

◆ event_adapt()

Adaptivity

◆ event_advection_term()

Advection term

◆ event_default_display()

◆ event_defaults()

Initial conditions

◆ event_init()

◆ event_projection()

◆ event_stability()

◆ event_viscous_term()

◆ prediction()

Predicted face velocity field

◆ properties()

◆ set_dtmax()

Time integration

◆ tracer_advection()

◆ tracer_diffusion()

◆ vof()

Variable Documentation

◆ a

◆ alpha

◆ dtmax

◆ g

◆ mgp

◆ mgpf

◆ mgu

◆ mu

◆ p

Incompressible Navier–Stokes solver (centered formulation)

◆ pf

◆ rho

◆ spread

◆ stokes

◆ uf

An advection solver