/** # 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](/src/references.bib#popinet2009) (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](http://www.ida.upmc.fr/~zaleski/paris/Vofi-1.0.tar.gz) library for accurate initialisation of interfacial shapes [Bna et al, 2015](/src/references.bib#bna2015). */ #include #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); } boundary ({c}); } /** 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() { 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); #endif /** Finally 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). */ fflush (stderr); system ("sort -k1,2 -n log > log.s; mv -f log.s log"); } /** 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](/src/references.bib#popinet2009). 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\$. ~~~gnuplot Relative curvature error as a function of resolution 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' ~~~ Similar results are obtained in three dimensions, but with a much wider intermediate range of low-order convergence for the max-norm. ~~~gnuplot Relative curvature error as a function of resolution in three dimensions 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' ~~~ */