src/test/neumann.c

Poisson equation on complex domains

We reproduce the test cases initially proposed by Johansen and Collela, 1998, Problem 1, p. 14, with Dirichlet boundary conditions and Problem 3, p. 19, with Neumann boundary conditions.

#include "embed.h"
#include "poisson.h"
#include "view.h"

The exact solution is given by ϕ(r,θ)=r4cos3θ

static double exact (double x, double y) {
  double θ = atan2 (y, x), r2 = x*x + y*y;
  return sq(r2)*cos (3.*θ);
}

double exact_gradient (Point point) {
  coord n = facet_normal (point, cs, fs);
  double nn = sqrt (sq(n.x) + sq(n.y));
  n.x /= nn, n.y /= nn;
  double θ = atan2(y,x), r = sqrt(x*x + y*y);
  double dsdtheta = - 3.*cube(r)*sin (3.*θ);
  double dsdr = 4.*cube(r)*cos (3.*θ);
  return (n.x*(dsdr*cos(θ) - dsdtheta*sin(θ)) +
	  n.y*(dsdr*sin(θ) + dsdtheta*cos(θ)));  
}

We also need a function for homogeneous Dirichlet condition.

static double dirichlet_homogeneous_bc (Point point, Point neighbor, scalar s) {
  return 0.;
}

int main()
{
  for (N = 32; N <= 512; N *= 2) {
    origin (-L0/2., -L0/2.);
    init_grid (N);

The shape of the domain is given by Ω=(r,θ):r0.30+0.15cos6θ for Problem 1 and Ω=(r,θ):r0.25+0.05cos6θ for Problem 2.

    vertex scalar φ[];
    foreach_vertex() {
      double θ = atan2(y, x), r = sqrt (sq(x) + sq(y));
#if DIRICHLET
      φ[] = 0.30 + 0.15*cos(6.*θ) - r;
#else
      φ[] = r - (0.25 + 0.05*cos(6.*θ));
#endif
    }
    fractions (φ, cs, fs);  
#if TREE
    cs.refine = cs.prolongation = fraction_refine;
#endif
    boundary ({cs,fs});
    restriction ({cs,fs});

    cm = cs;
    fm = fs;

Conditions on the box boundaries are set (only relevant for Problem 3).

    scalar a[], b[];
    
    a[left]   = exact (x - Δ/2., y);
    a.boundary_homogeneous[left] = dirichlet_homogeneous_bc;
    a[right]  = exact (x + Δ/2., y);
    a.boundary_homogeneous[right] = dirichlet_homogeneous_bc;
    a[top]    = exact (x, y + Δ/2.);
    a.boundary_homogeneous[top] = dirichlet_homogeneous_bc;
    a[bottom] = exact (x, y - Δ/2.);
    a.boundary_homogeneous[bottom] = dirichlet_homogeneous_bc;

The boundary conditions on the embedded boundary are Dirichlet and Neumann for Problem 1 and 3, respectively.

#if DIRICHLET
    a[embed] = dirichlet_embed (exact (x,y));
#else
    a[embed] = neumann_embed (exact_gradient (point));
#endif

The right-hand-side Δϕ=7r2cos3θ 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) : nodata;
      
      double xc = x, yc = y;
      if (cs[] > 0. && cs[] < 1.) {
	coord n = facet_normal (point, cs, fs), p;
	double α = plane_alpha (cs[], n);
	line_center (n, α, cs[], &p);
	xc += p.x*Δ, yc += p.y*Δ;
      }
      //      fprintf (stderr, "xc %g %g\n", xc, yc);
      double θ = atan2(yc, xc), r2 = sq(xc) + sq(yc);
      b[] = 7.*r2*cos (3.*θ)*cs[];
    }
    boundary ({a,b});

#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);

    // FIXME: need to set minlevel to 4
    timer t = timer_start();
    mgstats s = poisson (a, b, α = fs,
			 embed_flux =
			 a.boundary[embed] != symmetry ? embed_flux : NULL,
			 tolerance = 1e-6, minlevel = 4);
    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);
	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 is displayed using bview.

    view (fov = 18);
    squares ("a", spread = -1);
    draw_vof ("cs", "fs");
    save ("a.png");

    clear();
    squares ("e", spread = -1);
    draw_vof ("cs", "fs");
    save ("e.png");
#endif
    
    dump ("dump");
  }
}

Results

Problem 1

Solution on the finest mesh

Solution on the finest mesh

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

Problem 3

Solution on the finest mesh

Solution on the finest mesh

Error on the finest mesh

Error on the finest mesh

For Neumann boundary conditions, the residual also converges at first order in partial cells.

Maximum residual convergence

Maximum residual convergence

But this now leads to second-order overall convergence.

Maximum error convergence

Maximum error convergence

References

[johansen1998]

Hans Johansen and Phillip Colella. A cartesian grid embedded boundary method for poisson’s equation on irregular domains. Journal of Computational Physics, 147(1):60-85, 1998. [ .pdf ]