src/solve.h
Helper macro to invert (linear) spatial operators
The macro 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 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.
#include "poisson.h"Implementation
Note that the macro below is a slightly simplified version of the mg_solve() and mg_cycle() functions where more documentation can be found.
macro
mgstats solve (scalar a, double func, double rhs,
int nrelax = 4,
int minlevel = 0,
double tolerance = TOLERANCE)
{{
mgstats _s = (mgstats){0};
scalar _res[], _da[];
scalar_clone (_da, a);
for (int b = 0; b < nboundary; b++)
_da.boundary[b] = _da.boundary_homogeneous[b];
_s.nrelax = nrelax;
double _resb;
{
double maxres = 0.;
foreach (reduction(max:maxres)) {
_res[] = rhs - func;
if (fabs (_res[]) > maxres)
maxres = fabs (_res[]);
}
_resb = _s.resb = _s.resa = maxres;
}
for (_s.i = 0; _s.i < NITERMAX && (_s.i < NITERMIN || _s.resa > tolerance); _s.i++) {
{
restriction ({_res});
int _maxlevel = grid->maxdepth;
int _minlevel = min (minlevel, _maxlevel);
for (int l = _minlevel; l <= _maxlevel; l++) {
if (l == _minlevel)
foreach_level_or_leaf (l)
foreach_blockf (_da)
_da[] = 0.;
else
foreach_level (l)
foreach_blockf (_da)
_da[] = bilinear (point, _da);
boundary_level ({_da}, l);
for (int i = 0; i < _s.nrelax; i++) {
scalar a = _da;
foreach_level_or_leaf (l) {
a[] = 0.;
double _n = _res[] - func, _d;
diagonalize(a)
_d = func;
a[] = _n/_d;
}
boundary_level ({_da}, l);
}
}
foreach()
foreach_blockf (a)
a[] += _da[];
}
{
double maxres = 0.;
foreach (reduction(max:maxres)) {
_res[] = rhs - func;
if (fabs (_res[]) > maxres)
maxres = fabs (_res[]);
}
_s.resa = maxres;
}
if (_s.resa > tolerance) {
if (_resb/_s.resa < 1.2 && _s.nrelax < 100)
_s.nrelax++;
else if (_resb/_s.resa > 10 && _s.nrelax > 2)
_s.nrelax--;
}
_resb = _s.resa;
}
_s.minlevel = minlevel;
if (_s.resa > tolerance)
fprintf (stderr,
"src/solve.h:%d: warning: convergence for %s not reached after %d iterations\n"
" res: %g nrelax: %d\n", LINENO, a.name,
_s.i, _s.resa, _s.nrelax),
fflush (ferr);
return _s;
}}