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 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.

    vertex scalar φ[];
    foreach_vertex()
      φ[] = sq(0.392) - sq(x) - sq(y) - sq(z);
    boundary ({φ});
    fractions (φ, cs, fs);
#if TREE
    cs.refine = cs.prolongation = fraction_refine;
#endif
    boundary ({cs,fs});
    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 Δa=14a 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 α = plane_alpha (cs[], n);
	plane_center (n, α, cs[], &p);
	xc += p.x*Δ, yc += p.y*Δ, zc += p.z*Δ;
      }
      // fprintf (stderr, "xc %g %g %g\n", xc, yc, zc);
      b[] = - 14.*exact (xc, yc, zc)*cs[];
    }
    boundary ({a,b});

#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 α = plane_alpha (cs[], n);
	double length = line_length_center (n, α, &p);
	x += p.x*Δ, y += p.y*Δ;
	double θ = atan2(y,x), r = sqrt(x*x + y*y);
	
	double dsdtheta = - 3.*cube(r)*sin (3.*θ);
	double dsdr = 4.*cube(r)*cos (3.*θ);
	double nn = sqrt (sq(n.x) + sq(n.y));
	n.x /= nn, n.y /= nn;
	double dsdn = (n.x*(dsdr*cos(θ) - dsdtheta*sin(θ)) +
		       n.y*(dsdr*sin(θ) + dsdtheta*cos(θ)));

	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.α = fs;
    p.λ = 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, α = 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");
  }
}

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.

Maximum residual convergence

Maximum residual convergence

This leads to third-order overall convergence.

Error convergence

Error convergence

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).

Maximum residual convergence

Maximum residual convergence

But this now leads to second-order overall convergence.

Maximum error convergence

Maximum error convergence

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 ]