src/test/interp.c
Interpolation on halos
r(x,y)=sqrt(x*x+y*y)
set xrange [0:0.7]
set xlabel 'Radial coordinate'
set ylabel 'Error on halos'
plot '< awk "(\$3==10){print \$0}" out' u (r($1,$2)):($6) t 'level 10', \
'< awk "(\$3==9){print \$0}" out' u (r($1,$2)):($6) t 'level 9', \
'< awk "(\$3==8){print \$0}" out' u (r($1,$2)):($6) t 'level 8', \
'< awk "(\$3==7){print \$0}" out' u (r($1,$2)):($6) t 'level 7', \
'< awk "(\$3==6){print \$0}" out' u (r($1,$2)):($6) t 'level 6', \
'< awk "(\$3==5){print \$0}" out' u (r($1,$2)):($6) t 'level 5', \
'< awk "(\$3==4){print \$0}" out' u (r($1,$2)):($6) t 'level 4'
(script)
scalar h[];
int main (int argc, char ** argv)
{
int n = 2048;
init_grid (n);
double R0 = 0.1;
foreach()
h[] = exp(-(x*x + y*y)/(R0*R0));
/* initial coarsening (see halo.c) */
double tolerance = 1e-4;
while (adapt_wavelet ({h}, &tolerance, 11).nc);
#if 0
// we reinitialise h just to make sure that trash() does its job
trash ({h});
foreach()
h[] = exp(-(x*x + y*y)/(R0*R0));
#endif
double max = 0.;
for (int l = 0; l < depth(); l++)
foreach_halo (prolongation, l)
foreach_child() {
double e = exp(-(x*x+y*y)/(R0*R0)) - h[];
printf ("%g %g %d %d %g %g\n", x, y, level, cell.neighbors, h[], e);
if (fabs(e) > max)
max = fabs(e);
}
fprintf (stderr, "maximum error on halos: %g\n", max);
return (max > tolerance);
}