sandbox/Antoonvh/integrator.h
Integrate a field
#include "poisson.h"
double THETA_ANGLE = 0.0;
static double residual_int (scalar * al, scalar * bl, scalar * resl, void * data) {
double maxres = 0.;
double st = sin(THETA_ANGLE);
double ct = cos(THETA_ANGLE);
int ix = ct > 0 ? 0 : -1;
int iy = st > 0 ? 0 : -1;
scalar a = al[0], b = bl[0], res = resl[0];
foreach(reduction(max:maxres)) {
double ax = (a[1 + ix] - a[ix])/Delta;
double ay = (a[0, 1 + iy] - a[0, iy])/Delta;
res[] = -b[] - ct*ax - st*ay ;
if (fabs (res[]) > maxres)
maxres = fabs (res[]);
}
return maxres;
}
with a relevant relaxation-function formulation.
static void relax_int (scalar * al, scalar * bl, int l, void * data) {
scalar a = al[0];
scalar b = bl[0];
double st = sin(THETA_ANGLE);
double ct = cos(THETA_ANGLE);
int ix = ct > 0 ? 1 : -1;
int iy = st > 0 ? 1 : -1;
st = fabs(st);
ct = fabs(ct);
The scheme appears to benefit from a few relaxations iterations.
for (int j = 0; j < 2; j++){
foreach_level_or_leaf(l)
a[] = (-b[]*Delta + ct*a[ix] + st*a[0, iy])/(ct + st);
//a[] = (a[] + 2*(-b[]*Delta + ct*a[ix] + st*a[0, iy])/(ct + st))/3;
}
}
User interface
Let us hope this procedure converges.
struct Integrate_dx {
scalar a, b;
(const) face vector alpha;
(const) scalar lambda;
double tolerance;
int nrelax, minlevel;
scalar * res;
};
mgstats integrate_dx (struct Integrate_dx p) {
scalar a = p.a, b = p.b;
return mg_solve({a}, {b}, residual_int, relax_int,
&p, p.nrelax, p.res, minlevel = 1);
}