sandbox/Antoonvh/interface.c

The area of an interface

Basilisk provides and extensive toolbox to do interface reconstruction. On this page we check the convergence properties for determining the surface area of an iso-surface countour of an analytical function.

#include "grid/octree.h"
#include "fractions.h"
#include "utils.h"

More specifially, we try to numerically determine the surface area of a sphere with radius $R=0.5 $ a.u., such that the corresponding interface area should be A=\pi\ \text{ a.u.}^2. We use a function that varies quadtratically in space and one that varies linearly with the radius.

# define func (sq(x-xo) + sq(y-yo) + sq(z-zo) - sq(R)) 
# define func2 (pow(sq(x - xo) + sq(y - yo) + sq(z - zo), 0.5) - R)

double xo = 0.1, yo = M_PI/10.0;
double zo=-1./2.44, R = 0.5;
scalar f[], f2[];
int main(){
  L0 = 5.;
  X0 = Y0 = Z0 = -L0/2;
  init_grid(1 << 4);
  int i = 0;
  static FILE * fp1 = fopen("interface.dat", "w");
  for (i=1; i<8; i++){
    double s = 1./pow(2., (double)i); 
    refine(fabs(func) < (s) && level < i + 4);
    fraction(f, func);
    fraction(f2, func2);
    int cells = 0;
    foreach()
      if (f[] > 0. && f[] < 1.)
	cells++;
    fprintf(fp1, "%d\t%d\t%g\t%g\n",i, cells, fabs(interface_area(f)-M_PI), fabs(interface_area(f2)-M_PI));
    static FILE * fp = popen("gfsview-batch3D interface.interface.gfv | ppm2gif --delay 200 > surface.gif","w");
    output_gfs(fp);
    fprintf(fp, "Save stdout { format = PPM width = 600 height = 600}\n");
  }
  fclose(fp1);
}

First we visually inspect the effect of the increasing spatial resolution on the corresoponding fidelity in the represenation of the spherical surface.

Visual inspectation of the spatial convergence. The surface appears to be better resolved with each refinement step.

Next we check if the number of interfactial cells that we have used does indeed scale as expected.

set xr [0.5:7.5]
set logscale y
set ylabel 'number of interfacial cells' 
set xlabel 'refinement iteration'
set key top left
plot   50*2**(2*x) with lines title '2D scaling' ,\
      "interface.dat" using 1:2 title 'used interfacial cells'
    

2D scaling for a 2D problem. (script)

Now we check the error in the computed interfacial area:

set ylabel 'Error [a.u.^2]'
set key top right
plot "interface.dat" using 1:3 title 'Error in area f' ,\
     "interface.dat" using 1:4 title 'Error in area f2'

(script)

We see that it does not converge. This could be due to the fact that for the interface reconstruction from an analytical function, the fractions() function used (only) second-order-accurate interpolation to determine the ‘edge’ fraction of a cell’s edge. Thus when we integrate over O(n^2) cells, the second-order accuracy per cell results in O(n^2) \times O(n^{-2}) = O(N^0), i.e zeroth order global accuracy!