src/test/curvature.c
Curvature of a circular/spherical interface
This test evaluates the accuracy of the generalised height-function curvature calculation. It is similar to the case presented in Popinet, 2009 (Figure 5). The curvatures of circles/spheres with a randomised position and varying radii are computed and statistics on the error are gathered and displayed on the graph below.
We use the Vofi library for accurate initialisation of interfacial shapes Bna et al, 2015.
#include <vofi.h>
#pragma autolink -L$HOME/local/lib -lvofi
#include "fractions.h"
#include "curvature.h"We pass arguments to Vofi through these global variables.
static double xc, yc, zc, Rc;
static double sphere (creal p[dimension])
{
#if dimension == 3
return sq(p[0] - xc) + sq(p[1] - yc) + sq(p[2] - zc) - sq(Rc);
#else
return sq(p[0] - xc) + sq(p[1] - yc) - sq(Rc);
#endif
}
static void vofi (scalar c, int levelmax)
{
double fh = Get_fh (sphere, NULL, 1./(1 << levelmax), dimension, 0);
foreach() {
creal p[3] = {x - Delta/2., y - Delta/2., z - Delta/2.};
c[] = Get_cc (sphere, p, Delta, fh, dimension);
}
}The function below is called with different arguments for “coarse” and “fine” interface resolution in order to minimize the computation times. The number of randomised positions is nr, the radius of the circle is R, the maximum refinement level to consider is levelmax, and the statistics for each level are stored in the array n, while the statistics on which method is used are stored in sc.
void sample_circles (int nr, double R, int levelmax, norm * n, cstats * sc)
{
while (nr--) {We refine the grid down to levelmax but only around the interface.
scalar c[], kappa[];
c.refine = c.prolongation = fraction_refine;
init_grid (1 << 5);
xc = noise()/8., yc = noise()/8., zc = noise()/8., Rc = R;
vofi (c, 5);
for (int l = 6; l <= levelmax; l++) {
refine (c[] > 0. && c[] < 1. && level < l);
vofi (c, l);
}We then successively coarsen this fine initial grid to compute the curvature on coarser and coarser grids (thus saving on the expensive initial condition).
restriction ({c});
for (int l = levelmax; l >= 3; l--) {
unrefine (level >= l || c[] <= 0. || c[] >= 1.);
cstats s = curvature (c, kappa, sigma = 2.);We store statistics on the methods used for curvature computation…
sc[l].h += s.h; sc[l].f += s.f; sc[l].a += s.a; sc[l].c += s.c;
foreach()
if (c[] > 0. && c[] < 1.) {…and error statistics (for a given level of refinement l).
double e = fabs(kappa[]/2. - (dimension - 1)/R)*R/(dimension - 1);
n[l].volume += dv();
n[l].avg += dv()*e;
n[l].rms += dv()*e*e;
if (e > n[l].max)
n[l].max = e;
}
}
}
}
int main()
{
size (1. [0]); // dimensionless
origin (-0.5, -0.5, -0.5);
init_grid (N);We try a wide enough range of radii.
for (double R = 0.1; R <= 0.2; R *= 1.2) {We initialize the arrays required to store the statistics for each level of refinement.
int levelmax = 8;
norm n[levelmax + 1];
cstats sc[levelmax + 1];
for (int i = 0; i <= levelmax; i++) {
n[i].volume = n[i].avg = n[i].rms = n[i].max = 0;
sc[i].h = sc[i].f = sc[i].a = sc[i].c = 0.;
}We can limit randomisation for the higher resolutions (since we expect less “special cases” on fine meshes). We thus limit the total runtime by sampling many (100) locations on coarse meshes but only few (1) location on the finest mesh.
#if dimension == 2
sample_circles (1000, R, 4, n, sc);
sample_circles (100, R, 6, n, sc);
sample_circles (10, R, levelmax, n, sc);
#else // dimension == 3
sample_circles (100, R, 4, n, sc);
sample_circles (10, R, 6, n, sc);
sample_circles (1, R, levelmax, n, sc);
#endifFinally we output the statistics for this particular radius and for each level of refinement.
for (int l = levelmax; l >= 3; l--)
if (n[l].volume) {
n[l].avg /= n[l].volume;
n[l].rms = sqrt(n[l].rms/n[l].volume);
double t = sc[l].h + sc[l].f + sc[l].a + sc[l].c;
fprintf (stderr, "%g %g %g %g %g %g %g %g\n",
2.*R*(1 << l),
n[l].avg, n[l].rms, n[l].max,
sc[l].h/t, sc[l].f/t, sc[l].a/t, sc[l].c/t);
}
}At the end of the run, we sort the data by increasing order of diameter (in grid points).
The results are summarised in the figure below. There are two sets of points: error norms (bottom three curves) and percentages of curvatures computed with each method (top four curves). As expected the second-order convergence of the max and RMS norms is recovered for the pure HF method, for a diameter greater than about 15 grid points. The RMS norm shows consistent second-order convergence across the whole range of diameters. There is a significant degradation of the max norm between 7 and 15 grid points, corresponding with the introduction of the “averaging” and “(mixed) HF fit” methods. This degradation is significantly less pronounced for the method implemented in Popinet, 2009.
The grading of the methods used for curvature calculation follows what is expected:
- exclusively HF for D > 15,
- a combination of HF, nearest-neighbor average and “mixed HF fit” for 3 < D < 15
- and exclusively “centroids fit” for D < 3.
set logscale
set grid
set key top right
set xlabel 'Diameter (grid points)'
set ylabel 'Relative curvature error / percentage'
f(x)=(x > 0. ? 100.*x : 1e1000)
set yrange [:100]
plot 2./(x*x) t '2/x^{2}', 'log' u 1:4 w lp t 'Max', '' u 1:3 w lp t 'RMS', \
'../popinet.csv' u ($1*2):2 w lp t 'Popinet (2009)', \
'log' u 1:(f($5)) w lp t 'HF', '' u 1:(f($6)) w lp t 'HF fit', \
'' u 1:(f($7)) w lp t 'Average', '' u 1:(f($8)) w lp t 'Centroids'
Relative curvature error as a function of resolution (script)
Similar results are obtained in three dimensions, but with a much wider intermediate range of low-order convergence for the max-norm.
set key bottom left
plot 2./(x*x) t '2/x^{2}', '../curvature.3D/log' u 1:4 w lp t 'Max', \
'' u 1:3 w lp t 'RMS', '../curvature.3D/log' u 1:(f($5)) w lp t 'HF', \
'' u 1:(f($6)) w lp t 'HF fit', \
'' u 1:(f($7)) w lp t 'Average', '' u 1:(f($8)) w lp t 'Centroids'
Relative curvature error as a function of resolution in three dimensions (script)
