sandbox/Antoonvh/rf4t.c
Higher-order accurate face-value prolongation
Here we test the methods for 4th order face refinement and Solenoidal refinement.
#include "higher-order.h"
#include "poisson.h" //For projection()
#include "utils.h"
face vector fv[];
scalar s[];
double max_divergence (face vector fv);//Prototype
int main() {
L0 = 10;
X0 = Y0 = -L0/2;
foreach_dimension() {
periodic (left);
fv.x.prolongation = refine_face_4_x;
}
fv.x.refine = refine_face_solenoidal;
First, the prolongation method is tested:
set logscale xy 2
set xr [9:4000]
set yr [1e-12:1e-2]
set grid
set xlabel 'N'
set ylabel 'L_1 Error'
set size square
plot 'out' u 1:2, 1000*x**(-4)
for (int l = 4; l < 12; l++) {
init_grid (1 << l);
foreach_face() {
coord f = {y, -x};
fv.x[] = f.x*exp(-sq(x) - sq(y));
}
boundary ((scalar*){fv});
wavelet (fv.x, s);
fprintf (stdout, "%d %g\n", N, normf(s).avg);
}
Next, the Solenoidal refinement test is performed
Accuracy:
set logscale xy 2
set xr [9:4000]
set yr [1e-13:1e-2]
set grid
set size square
set xlabel 'N'
set ylabel 'L_1 Error'
plot 'log' u 1:2, 500000*x**(-5)
Divergence:
reset
set logscale x
set xr [9:4000]
set yr [-1e-11:1e-11]
set grid
set xlabel 'N'
set ylabel 'Max. absolute divergence'
set size square
plot 'log' u 1:4
NITERMIN = 25;
for (int l = 4; l < 12; l++) {
init_grid (1 << l);
foreach_face() {
coord f = {y, -x};
fv.x[] = f.x*exp(-sq(x) - sq(y));
}
project (fv, s);
double md = max_divergence (fv); //approx zero
double * fva = malloc (2*(sq(N))*sizeof(double));
long int j = 0;
foreach_face()
fva[j++] = fv.x[]; //Array of reference values.
unrefine (level == l - 1);
refine (level < l);
double err = 0;
j = 0;
foreach_face()
err += fabs(fv.x[] - fva[j++]);
fprintf (stderr, "%d %g %g %g\n",
N, err/sq(N), md, max_divergence (fv));
free (fva);
}
}
double max_divergence (face vector fv) {
double max = 0;
foreach() {
double d = 0;
foreach_dimension()
d += fv.x[1] - fv.x[];
d /= Delta;
if (fabs(d) > max)
max = fabs(d);
}
return max;
}