sandbox/Antoonvh/refineing2d.c

The accuracy of various refinement/prolongation techniques in 2D

Here we test the accuracy of various techniques used for Prolongation and/or refinement.

The algorithm

We determine the accuracy of a a few interpolation techniques by applying them to a 2D Gaussian distribution. We therefore use the quadtree toolbox, the statsf function and two additional attributes.

#include "grid/quadtree.h"
#include "utils.h"
#include "additional_attributes.h"

scalar c[],err[];
FILE * fp1;
int cells = 0;
stats s;
int maxlevel = 11;

int main(){
  fp1=fopen("data.dat","w");
  L0=10.;
  X0=Y0=-L0/2;

We start with initializing the Gaussian distribution on a fine-resolution grid corresponding to maxlevel.

  init_grid(1<<maxlevel);
  foreach()
    c[]=exp(-(sq(x)+sq(y)));

In this loop the analysis is done on succesively coarsened grid.

  for (int n=1;n<8;n++){
    cells=0;
    foreach()
      cells++;
    fprintf(ferr,"%d\t %d\n",n,cells);
    fprintf(fp1,"%d\t%g\t",n,L0/(pow(2.,(double)(maxlevel-n+1.))));

For each resolution, five different interpolation techniques are tested.

    for (int j=1;j<7;j++){
      if (j==1) c.prolongation=refine_injection;
      if (j==2) c.prolongation=refine_bilinear;
      if (j==3) c.prolongation=refine_linear;
      if (j==4) c.prolongation=refine_quadratic; //Third order?
      if (j==5) c.prolongation=refine_3x3_linear; // different from linear?
      if (j==6) c.prolongation=refine_biquadratic; //Third order?

The error is evaluated using the wavelet methodology.

      wavelet(c,err);
      foreach()
	err[]=fabs(err[]);
      s=statsf(err);
      fprintf(fp1,"%g\t",s.sum);
    }
    fprintf(fp1,"\n");

After the results are written to a file, the grid is coarsened such that the analysis can be repeated with a different \Delta.

    unrefine(level>=(maxlevel-n)); 
  }
}

Results

  set xr [0.004:0.5]
  set yr [0.000001:1]
  set logscale y
  set logscale x
  set xlabel '{/Symbol D}'
  set ylabel 'Sum of Error'
  set key right bottom box 1
  set size square
  plot   (2*x**1) lw 3 lc rgb 'purple' title '{/Symbol \265}{/Symbol D}^{1}' ,\
   (2*x**2) lw 3 lc rgb 'blue' title '{/Symbol \265}{/Symbol D}^{2}' ,\
    (2*x**3) lw 3 lc rgb 'red' title '{/Symbol \265}{/Symbol D}^{3}' ,\
  'data.dat' using 2:3 title 'Injection',\
  'data.dat' using 2:4 title 'Bilinear' ,\
  'data.dat' using 2:5 title 'Linear' ,\
  'data.dat' using 2:6 title 'quadratic' ,\
  'data.dat' using 2:7 title '3x3-Linear' ,\
  'data.dat' using 2:8 title 'biquadratic'
  

Results from various interpolation/upsample techniques (script)

Appearantly the Injection method is first-order accurate, Linear and Bilinear are both second-order accurate and the 3x3-linear formulation is third-order accurate. The quadratic interpolation technique, that was designed to fall in 3-rd order catogory, turns out to be only second-order accurate. So something went wrong there.

Also note that the Injection, Linear and 3x3-linear techniques are conservative for the first-order moments. This could be favourable, as it may cause the numerical formulation of refinement to enherit some of the properties from the physical system of your interest (e.g. the refinement of a volicity field).

Todo

• Find out why the quadratic method is not third-order accurate
• Implement a 3x3(x3)-linear analogue for 3D fields and test this if it is also third-order accurate. (DONE see elswere)