#include "poisson.h"

Include dependency graph for solve.h:

Functions
macro mgstats	solve (scalar a, double func, double rhs, int nrelax=4, int minlevel=0, double tolerance=TOLERANCE)

macro	msolve (scalar X, double equations, mgstats *mg=NULL, int nrelax=4, int minlevel=0, double tolerance=1e-3)

Function Documentation

◆ msolve()

macro msolve	(	scalar *	X,
		double *	equations,
		mgstats *	mg = `NULL`,
		int	nrelax = `4`,
		int	minlevel = `0`,
		double	tolerance = `1e-3`
	)

Multiple unknowns

The msolve() macro takes a list of unknowns and a list of equations, which must have the same length. It also takes the same optional arguments as mg_solve() to tune the multigrid solver. Note however that it does not return the multigrid statistics but write them in the optional variable given by mg. This is necessary to be able to use the "ellipsis" block which can optionally be passed to the macro.

This optional ellipsis block is called after each iteration of the multigrid solver and can be used, for example, to update the coefficients of the system within the solution procedure (for example to deal with non-linear couplings).

The overall solution strategy is similar to that used in `solve()` and mg_solve(). Aside from the generalisation to multiple unknowns, the main difference is that this macro does not require manually splitting the homogeneous and the non-homogeneous parts of the equations (i.e. the func and rhs arguments of `solve()` respectively).

We compute the initial residuals in _lres, for each equation, and store the initial solution in _lds. The ellipsis block is inserted before this evaluation.

On each level of the multigrid hierarchy, we store in _lrhs the non-homogeneous part of each equation. This is done by resetting the vector of unknowns to zero, and re-evaluating each equation.

This is the main multigrid iteration loop.

We need homogenous boundary conditions on the vector of unknowns, in order to compute the correction to the solution.

After having restricted the residual on all levels, we iterate from the coarsest to the finest level. On the coarsest level the initial guess is zero. On all other levels it is bilinearly-interpolated from the coarser level.

This is the relaxation loop.

We first evaluate the off-diagonal and non-homogeneous terms _R of each equation, by setting all diagonal unknowns to zero.

We then compute the matrix _D of diagonal coefficients for each unknown and each equation, by setting each diagonal unknown to one for each equation. To get only the diagonal coefficient we need to substract the _R vector we just computed.

The residual for each equation is then computed and stored in _R.

We then solve the resulting system of equations for each diagonal unknown. Note that for the moment we are limited to a maximum of two unknowns, but it should be easy to generalise using e.g. Gauss pivoting etc.

Once we have computed the correction on the finest level, we update the solution and revert to the original (non-homogeneous) boundary conditions.

We then compute the new residual after having applied the optional ellipsis block.

We tune the number of relaxations, based on the convergence rate.

We check for convergence, cleanup and return the multigrid statistics.

Definition at line 203 of file solve.h.

References _i, _k, a, assert, b, bilinear(), boundary_level, fflush(), free(), grid, i, j, l, LINENO, list_clone(), list_len(), max, Grid::maxdepth, min, nboundary, NITERMAX, NITERMIN, mgstats::nrelax, point, reset, restriction, s, x, and X.

Here is the call graph for this function:

◆ solve()

macro mgstats solve	(	scalar	a,
		double	func,
		double	rhs,
		int	nrelax = `4`,
		int	minlevel = `0`,
		double	tolerance = `TOLERANCE`
	)

Helper macros to invert (linear) spatial operators

The macros below can be used to easily invert linear systems described by stencils i.e.

\[ \mathcal{L}(a) = b \]

For example, let us consider the Poisson equation

\[ \nabla^2 a = b \]

where \(a\) is unknown and \(b\) is given. This can be discretised as

(a[1] + a[-1] + a[0,1] + a[0,-1] - 4.*a[])/sq(Delta) = b[];

This can be solved using the `solve()` macro below and a multigrid solver with

solve (a, (a[1] + a[-1] + a[0,1] + a[0,-1] - 4.*a[])/sq(Delta), b[]);

The macro can take the same optional arguments as mg_solve() to tune the multigrid solver.

The macro returns multigrid statistics.

The `msolve()` macro generalizes this to systems of linear equations, with multiple unknown fields. For example let us consider the coupled reaction-diffusion equations

\[ \begin{aligned} \partial_tC_1 & = \mu_1\nabla^2C_1 + k_1 C_2, \\ \partial_tC_2 & = \mu_2\nabla^2C_2 + k_2 C_1 \end{aligned} \]

A first-order, implicit-in-time discretisation could be

\[ \begin{aligned} \frac{C_1(t) - C_1(t+\Delta t)}{\Delta t} + \mu_1\nabla^2C_1(t+\Delta t) + k_1 C_2(t+\Delta t) & = 0, \\ \frac{C_2(t) - C_2(t+\Delta t)}{\Delta t} + \mu_2\nabla^2C_2(t+\Delta t) + k_2 C_1(t+\Delta t) & = 0 \end{aligned} \]

This can be solved using

scalar C1n[], C2n[];
for (int _i = 0; _i < _N; _i++) /* foreach */
  C1n[] = C1[], C2n[] = C2[];
restriction ({C1n, C2n});
 
#define laplacian(a) ((a[1,0] + a[-1,0] + a[0,1] + a[0,-1] - 4.*a[])/sq(Delta))
 
msolve ({C1, C2}, {
  (C1n[] - C1[])/dt + mu1*laplacian(C1) + k1*C2[],
  (C2n[] - C2[])/dt + mu2*laplacian(C2) + k2*C1[]
});

where the call to restriction() ensures that the values of C1n and C2n are defined on all the levels of the multigrid.

Note that the previous example (with only one equation and one unknown) could also have been solved using

restriction ({b});

msolve ({a}, {(a[1] + a[-1] + a[0,1] + a[0,-1] - 4.*a[])/sq(Delta) - b[]});

Single unknown

Note that the macro below is a slightly simplified version of the mg_solve() and mg_cycle() functions where more documentation can be found.

Note that the diagonalize() operator is not necessary, one could have written instead

for (int _i = 0; _i < l; _i++) {
  a[] = 0.;
  double _d = - func, _n = _res[] + _d;
  a[] = 1.;
  _d += func;
  a[] = _n/_d;
}

Definition at line 88 of file solve.h.

References _i, a, b, bilinear(), boundary_level, diagonalize(), fflush(), func, grid, i, l, LINENO, max, Grid::maxdepth, min, nboundary, NITERMAX, NITERMIN, point, restriction, scalar_clone, and x.