src/test/neumann3D.c
Poisson equation on complex domains in 3D
We reproduce the test cases initially proposed by Schwarz et al., 2006, section 4.1, with Dirichlet or Neumann boundary conditions.
#include "grid/multigrid3D.h"
#include "embed.h"
#include "poisson.h"
#include "view.h"
The exact solution is given by \displaystyle a(x,y,z) = \sin(x) \sin(2y) \sin(3z)
static double exact (double x, double y, double z) {
return sin(x)*sin(2.*y)*sin(3.*z);
}
double exact_gradient (Point point, double xc, double yc, double zc)
{
coord n = facet_normal (point, cs, fs);
normalize (&n);
return (n.x*cos(xc)*sin(2.*yc)*sin(3.*zc) +
n.y*2.*sin(xc)*cos(2.*yc)*sin(3.*zc) +
n.z*3.*sin(xc)*sin(2.*yc)*cos(3.*zc));
}
int main()
{
for (N = 32; N <= 256; N *= 2) {
size (1. [0]); // dimensionless space
origin (-L0/2., -L0/2., -L0/2.);
init_grid (N);
The domain is a sphere of radius 0.392 centered at the origin.
solid (cs, fs, sq(0.392) - sq(x) - sq(y) - sq(z));
#if TREE
cs.refine = cs.prolongation = fraction_refine;
#endif
restriction ({cs,fs});
cm = cs;
fm = fs;
scalar a[], b[];
The boundary conditions on the embedded boundary are Dirichlet and Neumann, respectively.
#if DIRICHLET
a[embed] = dirichlet (exact (x, y, z));
#else
a[embed] = neumann (exact_gradient (point, x, y, z));
#endif
We use “third-order” face flux interpolation.
a.third = true;
The right-hand-side \displaystyle \Delta a = -14 a is defined using the coordinates of the barycenter of the cut cell (xc,yc), which is calculated from the cell and surface fractions.
foreach() {
a[] = cs[] > 0. ? exact (x, y, z) : nodata;
double xc = x, yc = y, zc = z;
if (cs[] > 0. && cs[] < 1.) {
coord n = facet_normal (point, cs, fs), p;
double alpha = plane_alpha (cs[], n);
plane_center (n, alpha, cs[], &p);
xc += p.x*Delta, yc += p.y*Delta, zc += p.z*Delta;
}
// fprintf (stderr, "xc %g %g %g\n", xc, yc, zc);
b[] = - 14.*exact (xc, yc, zc)*cs[];
}
#if 0
foreach_face (z)
if (z == 0. && fs.z[] > 0. && fs.z[] < 1.) {
// fprintf (stderr, "fs %g %g %g %g\n", x, y, z, fs.z[]);
embed_face_gradient_z (point, a, 0);
}
exit (0);
#endif
#if 0
output_cells (stdout);
output_facets (cs, stdout, fs);
scalar e[];
foreach() {
if (cs[] > 0. && cs[] < 1.) {
scalar s = a;
coord n = facet_normal (point, cs, fs), p;
double alpha = plane_alpha (cs[], n);
double length = line_length_center (n, alpha, &p);
x += p.x*Delta, y += p.y*Delta;
double theta = atan2(y,x), r = sqrt(x*x + y*y);
double dsdtheta = - 3.*cube(r)*sin (3.*theta);
double dsdr = 4.*cube(r)*cos (3.*theta);
double nn = sqrt (sq(n.x) + sq(n.y));
n.x /= nn, n.y /= nn;
double dsdn = (n.x*(dsdr*cos(theta) - dsdtheta*sin(theta)) +
n.y*(dsdr*sin(theta) + dsdtheta*cos(theta)));
e[] = dsdn - dirichlet_gradient (point, s, cs, n, p, exact (x, y));
#if 1
fprintf (stderr, "g %g %g %g %g\n",
x, y, dsdn,
dirichlet_gradient (point, s, cs, n, p, exact (x, y)));
#endif
}
else
e[] = nodata;
}
norm n = normf (e);
fprintf (stderr, "%d %g %g\n",
N, n.rms, n.max);
#else
The Poisson equation is solved.
struct Poisson p;
p.alpha = fs;
p.lambda = zeroc;
p.embed_flux = embed_flux;
scalar res[];
double maxp = residual ({a}, {b}, {res}, &p), maxf = 0.;
foreach()
if (cs[] == 1. && fabs(res[]) > maxf)
maxf = fabs(res[]);
fprintf (stderr, "maxres %d %.3g %.3g\n", N, maxf, maxp);
timer t = timer_start();
mgstats s = poisson (a, b, alpha = fs,
flux = a.boundary[embed] != symmetry ? embed_flux : NULL,
tolerance = 1e-6);
double dt = timer_elapsed (t);
printf ("%d %g %d %d\n", N, dt, s.i, s.nrelax);
The total (e), partial cells (ep) and full cells (ef) errors fields and their norms are computed.
scalar e[], ep[], ef[];
foreach() {
if (cs[] == 0.)
ep[] = ef[] = e[] = nodata;
else {
e[] = a[] - exact (x, y, z);
ep[] = cs[] < 1. ? e[] : nodata;
ef[] = cs[] >= 1. ? e[] : nodata;
}
}
norm n = normf (e), np = normf (ep), nf = normf (ef);
fprintf (stderr, "%d %.3g %.3g %.3g %.3g %.3g %.3g %d %d\n",
N, n.avg, n.max, np.avg, np.max, nf.avg, nf.max, s.i, s.nrelax);
The solution error is displayed using bview.
view (fov = 16.1659, quat = {-0.270921,0.342698,0.106093,0.893256},
tx = -0.00535896, ty = 0.000132663, bg = {1,1,1},
width = 600, height = 600, samples = 4);
draw_vof("cs", "fs", color = "e");
save ("e.png");
#endif
// dump ("dump"); // too big
}
}
Results
Dirichlet boundary condition

Error on the finest mesh
For Dirichlet boundary conditions, the residual converges at first order in partial cells.
set xrange [*:*]
ftitle(a,b) = sprintf("%.3f/x^{%4.2f}", exp(a), -b)
f(x) = a + b*x
fit f(x) '< grep maxres ../dirichlet3D/log' u (log($2)):(log($3)) via a,b
f2(x) = a2 + b2*x
fit f2(x) '' u (log($2)):(log($4)) via a2,b2
set xlabel 'Resolution'
set logscale
set cbrange [1:2]
set xtics 16,2,1024
set grid ytics
set ytics format "% .0e"
set xrange [16:512]
set ylabel 'Maximum residual'
set yrange [1e-6:]
set key bottom left
plot '' u 2:3 t 'full cells', exp(f(log(x))) t ftitle(a,b), \
'' u 2:4 t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2)
Maximum residual convergence (script)
This leads to third-order overall convergence.
set xrange [*:*]
fit f(x) '../dirichlet3D/log' u (log($1)):(log($3)) via a,b
fit f2(x) '' u (log($1)):(log($4)) via a2,b2
f3(x) = a3 + b3*x
fit f3(x) '' u (log($1)):(log($6)) via a3,b3
set xrange [16:512]
set ylabel 'Error'
set yrange [*:*]
plot '' u 1:3 pt 6 t 'max all cells', exp(f(log(x))) t ftitle(a,b), \
'' u 1:4 t 'avg partial cells', exp(f2(log(x))) t ftitle(a2,b2), \
'' u 1:6 t 'avg full cells', exp(f3(log(x))) t ftitle(a3,b3)
Error convergence (script)
Neumann boundary condition

Error on the finest mesh
For Neumann boundary conditions, the residual converges at less than first order in partial cells (which may be worth looking into).
set xrange [*:*]
fit f(x) '< grep maxres log' u (log($2)):(log($3)) via a,b
fit f2(x) '' u (log($2)):(log($4)) via a2,b2
set ylabel 'Maximum residual'
set xrange [16:512]
set yrange [1e-6:]
set key bottom left
plot '' u 2:3 t 'full cells', exp(f(log(x))) t ftitle(a,b), \
'' u 2:4 t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2)
Maximum residual convergence (script)
But this now leads to second-order overall convergence.
set xrange [*:*]
fit f(x) 'log' u (log($1)):(log($7)) via a,b
fit f2(x) '' u (log($1)):(log($5)) via a2,b2
set ylabel 'Maximum error'
set xrange [16:512]
set yrange [*:*]
plot '' u 1:7 pt 6 t 'full cells', exp(f(log(x))) t ftitle(a,b), \
'' u 1:5 t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2)
Maximum error convergence (script)
References
[schwartz2006] |
Peter Schwartz, Michael Barad, Phillip Colella, and Terry Ligocki. A cartesian grid embedded boundary method for the heat equation and poisson’s equation in three dimensions. Journal of Computational Physics, 211(2):531–550, 2006. [ .pdf ] |