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

Data Structures
struct	mgstats
	Information about the convergence of the solver is returned in a structure. More...

struct	Poisson

Functions
trace void	mg_cycle (scalar a, scalar res, scalar da, void(relax)(scalar da, scalar res, int depth, void data), void data, int nrelax, int minlevel, int maxlevel)

trace mgstats	mg_solve (scalar a, scalar b, double(residual)(scalar a, scalar b, scalar res, void data), void(relax)(scalar da, scalar res, int depth, void data), void data=NULL, int nrelax=4, scalar *res=NULL, int minlevel=0, double tolerance=TOLERANCE)
	The user needs to provide a function which computes the residual field (and returns its maximum) as well as the relaxation function.

static void	relax (scalar al, scalar bl, int l, void *data)
	We can now write the relaxation function.

static double	residual (scalar al, scalar bl, scalar resl, void data)
	The equivalent residual function is obtained in a similar way in the case of a Cartesian grid, however the case of the tree mesh requires more careful consideration...

trace mgstats	poisson (scalar a, scalar b, const vector alpha={{-1}}, const scalar lambda={-1}, double tolerance=0., int nrelax=4, int minlevel=0, scalar res=NULL, double(flux)(Point, scalar, vector, double *)=NULL)

trace mgstats	project (vector uf, scalar p, const vector alpha=unityf, double dt=1., int nrelax=4)

Variables
int	NITERMAX = 100

int	NITERMIN = 1

double	TOLERANCE = 1e-3 [*]

Function Documentation

◆ mg_cycle()

trace void mg_cycle	(	scalar *	a,
		scalar *	res,
		scalar *	da,
		void()(scalar da, scalar res, int depth, void data)	relax,
		void *	data,
		int	nrelax,
		int	minlevel,
		int	maxlevel
	)

Multigrid Poisson–Helmholtz solvers

We want to solve Poisson–Helmholtz equations of the general form

\[ L(a) = \nabla\cdot (\alpha\nabla a) + \lambda a = b \]

This can be done efficiently using a multigrid solver.

An important aspect of Poisson–Helmholtz equations is that the operator \(L()\) is linear. This property can be used to build better estimates of a solution by successive corrections to an initial guess. If we define an approximate solution \(\tilde{a}\) as

\[ \tilde{a} + da = a \]

where \(a\) is the exact (unknown) solution, using the linearity of the operator we find that \(da\) verifies

\[ L(da) = b - L(\tilde{a}) \]

where the right-hand-side is often called the residual of the approximate solution \(\tilde{a}\).

Multigrid cycle

Here we implement the multigrid cycle proper. Given an initial guess a*, a residual res, a correction field da and a relaxation function relax, we will provide an improved guess at the end of the cycle.

We first define the residual on all levels.

We then proceed from the coarsest grid (minlevel) down to the finest grid.

On the coarsest grid, we take zero as initial guess.

On all other grids, we take as initial guess the approximate solution on the coarser grid bilinearly interpolated onto the current grid.

We then apply homogeneous boundary conditions and do several iterations of the relaxation function to refine the initial guess.

And finally we apply the resulting correction to a.

Definition at line 35 of file poisson.h.

References _i, a, bilinear(), boundary_level, data, i, l, min, point, relax(), restriction, s, and x.

Referenced by mg_solve().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ mg_solve()

trace mgstats mg_solve	(	scalar *	a,
		scalar *	b,
		double()(scalar a, scalar b, scalar res, void *data)	residual,
		void()(scalar da, scalar res, int depth, void data)	relax,
		void *	data = `NULL`,
		int	nrelax = `4`,
		scalar *	res = `NULL`,
		int	minlevel = `0`,
		double	tolerance = `TOLERANCE`
	)

The user needs to provide a function which computes the residual field (and returns its maximum) as well as the relaxation function.

The user-defined pointer data can be used to pass arguments to these functions. The optional number of relaxations is nrelax and res is an optional list of fields used to store the residuals. The minimum level of the hierarchy can be set (default is zero i.e. the root cell).

We allocate a new correction and residual field for each of the scalars in a.

The boundary conditions for the correction fields are the homogeneous* equivalent of the boundary conditions applied to a*.

We initialise the structure storing convergence statistics.

Here we compute the initial residual field and its maximum.

We then iterate until convergence or until NITERMAX is reached. Note also that we force the solver to apply at least one cycle, even if the initial residual is lower than TOLERANCE.

We tune the number of relaxations so that the residual is reduced by between 2 and 20 for each cycle. This is particularly useful for stiff systems which may require a larger number of relaxations on the finest grid.

If we have not satisfied the tolerance, we warn the user.

We deallocate the residual and correction fields and free the lists.

Definition at line 130 of file poisson.h.

References a, b, data, fflush(), free(), grid, scalar::i, LINENO, list_clone(), Grid::maxdepth, mg_cycle(), nboundary, NITERMAX, NITERMIN, p, relax(), s, sum, v, and x.

Referenced by axistream(), event_pressure(), event_tracer_diffusion(), poisson_thermal(), update_green_naghdi(), and viscosity().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ poisson()

trace mgstats poisson	(	scalar	a,
		scalar	b,
		const vector	alpha = `{{-1}}`,
		const scalar	lambda = `{-1}`,
		double	tolerance = `0.`,
		int	nrelax = `4`,
		int	minlevel = `0`,
		scalar *	res = `NULL`,
		double()(Point, scalar, vector, double )	flux = `NULL`
	)

User interface

Finally we provide a generic user interface for a Poisson–Helmholtz equation of the form

\[ \nabla\cdot (\alpha\nabla a) + \lambda a = b \]

If \(\alpha\) or \(\lambda\) are not set, we replace them with constant unity vector (resp. zero scalar) fields. Note that the user is free to provide \(\alpha\) and \(\beta\) as constant fields.

We need \(\alpha\) and \(\lambda\) on all levels of the grid.

If tolerance is set it supersedes the default of the multigrid solver.

We restore the default.

Definition at line 411 of file poisson.h.

◆ project()

trace mgstats project	(	vector	uf,
		scalar	p,
		const vector	alpha = `unityf`,
		double	dt = `1.`,
		int	nrelax = `4`
	)

Projection of a velocity field

The function below "projects" the velocity field u onto the space of divergence-free velocity fields i.e.

\[ \mathbf{u}_f^{n+1} \leftarrow \mathbf{u}_f - \Delta t\alpha\nabla p \]

so that

\[ \nabla\cdot\mathbf{u}_f^{n+1} = 0 \]

This gives the Poisson equation for the pressure

\[ \nabla\cdot(\alpha\nabla p) = \frac{\nabla\cdot\mathbf{u}_f}{\Delta t} \]

We allocate a local scalar field and compute the divergence of \(\mathbf{u}_f\). The divergence is scaled by dt so that the pressure has the correct dimension.

We solve the Poisson problem. The tolerance (set with TOLERANCE) is the maximum relative change in volume of a cell (due to the divergence of the flow) during one timestep i.e. the non-dimensional quantity

\[ |\nabla\cdot\mathbf{u}_f|\Delta t \]

Given the scaling of the divergence above, this gives

And compute \(\mathbf{u}_f^{n+1}\) using \(\mathbf{u}_f\) and \(p\).

Definition at line 483 of file poisson.h.

References _i, alpha, dimension, dt, face_gradient_x, mgp, Poisson::nrelax, p, poisson, sq(), TOLERANCE, Poisson::tolerance, uf, vector::x, and x.

Referenced by event_advection_term(), event_end_timestep(), and event_projection().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ relax()

static void relax	(	scalar *	al,
		scalar *	bl,
		int	l,
		void *	data
	)

static

We can now write the relaxation function.

We first recover the extra parameters from the data pointer.

We use either Jacobi (under)relaxation, Gauss-Seidel or we directly reuse values as soon as they are updated. For Jacobi, we need to allocate space for the new field c. Jacobi is useful mostly as it gives results which are independent of the order in which the cells are traversed. This is not the case for the simple traversal, which means for example that results will depend on whether a tree or a multigrid is used (because cells will be traversed in a different order). The same comment applies to OpenMP or MPI parallelism. In practice however Jacobi convergence tends to be slower than simple reuse.

On GPUs, we use red/black Gauss-Seidel relaxation, which requires two loops (for odd/even indices). Note also that, unlike the other option, red/black relaxation should be deterministic.

We use the face values of \(\alpha\) to weight the gradients of the 5-points Laplacian operator. We get the relaxation function.

For weighted Jacobi we under-relax with a weight of 2/3.

Definition at line 272 of file poisson.h.

References _i, a, alpha, b, c, d, data, dimension, Point::i, l, lambda, level, n, p, point, sq(), vector::x, and x.

Referenced by mg_cycle(), and mg_solve().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ residual()

static double residual	(	scalar *	al,
		scalar *	bl,
		scalar *	resl,
		void *	data
	)

static

The equivalent residual function is obtained in a similar way in the case of a Cartesian grid, however the case of the tree mesh requires more careful consideration...

Definition at line 356 of file poisson.h.

References _i, a, alpha, b, c, data, dimension, face_gradient_x, g, lambda, max, p, point, Poisson::res, vector::x, and x.

Variable Documentation

◆ NITERMAX

int NITERMAX = 100

Multigrid solver

The multigrid solver itself uses successive calls to the multigrid cycle to refine an initial guess until a specified tolerance is reached.

The maximum number of iterations is controlled by NITERMAX and the tolerance by TOLERANCE with the default values below.

Definition at line 106 of file poisson.h.

Referenced by mg_solve(), msolve(), and solve().

◆ NITERMIN

int NITERMIN = 1

Definition at line 106 of file poisson.h.

Referenced by mg_solve(), msolve(), and solve().

◆ TOLERANCE

double TOLERANCE = 1e-3 [*]

Definition at line 107 of file poisson.h.

Referenced by event_pressure(), event_update_eta(), poisson_thermal(), project(), and viscosity().

Data Structures

Functions

Variables

Function Documentation

◆ mg_cycle()

Multigrid Poisson–Helmholtz solvers

Multigrid cycle

◆ mg_solve()

◆ poisson()

User interface

◆ project()

Projection of a velocity field

◆ relax()

◆ residual()

Variable Documentation

◆ NITERMAX

Multigrid solver

◆ NITERMIN

◆ TOLERANCE