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) {
    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,
			 embed_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

    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)

    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)

    Error convergence (script)

    Neumann boundary condition

    Error on the finest mesh

    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)

    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)

    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 ]