# 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[]);
}
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

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)

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)

### Neumann boundary condition

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)

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)

## 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 ]