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 \displaystyle \phi(r,\theta) = r^4\cos 3\theta

    static double exact (double x, double y) {
      double theta = atan2 (y, x), r2 = x*x + y*y;
      return sq(r2)*cos (3.*theta);
    }
    
    double exact_gradient (Point point, double theta, double r)
    {
      coord n = facet_normal (point, cs, fs);
      normalize (&n);
      double dsdtheta = - 3.*cube(r)*sin (3.*theta);
      double dsdr = 4.*cube(r)*cos (3.*theta);
      return (n.x*(dsdr*cos(theta) - dsdtheta*sin(theta)) +
    	  n.y*(dsdr*sin(theta) + dsdtheta*cos(theta)));  
    }

    We also need a function for homogeneous Dirichlet condition.

    static
    double dirichlet_homogeneous_bc (Point point, Point neighbor,
    				 scalar s, void * data) {
      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 \displaystyle \Omega = {(r,\theta): r \leq 0.30 + 0.15\cos 6\theta} for Problem 1 and \displaystyle \Omega = {(r,\theta): r \geq 0.25 + 0.05\cos 6\theta} for Problem 2.

        vertex scalar phi[];
        foreach_vertex() {
          double theta = atan2(y, x), r = sqrt (sq(x) + sq(y));
    #if DIRICHLET
          phi[] = 0.30 + 0.15*cos(6.*theta) - r;
    #else
          phi[] = r - (0.25 + 0.05*cos(6.*theta));
    #endif
        }
        boundary ({phi});
        fractions (phi, 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 - Delta/2., y);
        a.boundary_homogeneous[left] = dirichlet_homogeneous_bc;
        a[right]  = exact (x + Delta/2., y);
        a.boundary_homogeneous[right] = dirichlet_homogeneous_bc;
        a[top]    = exact (x, y + Delta/2.);
        a.boundary_homogeneous[top] = dirichlet_homogeneous_bc;
        a[bottom] = exact (x, y - Delta/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 (exact (x,y));
    #else
        a[embed] = neumann (exact_gradient (point, atan2(y, x), sqrt(x*x + y*y)));
    #endif

    We use “third-order” face flux interpolation.

        a.third = true;

    The right-hand-side \displaystyle \Delta\phi = 7r^2\cos 3\theta 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 alpha = plane_alpha (cs[], n);
    	line_center (n, alpha, cs[], &p);
    	xc += p.x*Delta, yc += p.y*Delta;
          }
          //      fprintf (stderr, "xc %g %g\n", xc, yc);
          double theta = atan2(yc, xc), r2 = sq(xc) + sq(yc);
          b[] = 7.*r2*cos (3.*theta)*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 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);
    
        // FIXME: need to set minlevel to 4
        timer t = timer_start();
        mgstats s = poisson (a, b, alpha = fs, 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.

    set xrange [*:*]
    ftitle(a,b) = sprintf("%.3f/x^{%4.2f}", exp(a), -b)
    f(x) = a + b*x
    fit f(x) '< grep maxres ../dirichlet/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:1024]
    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), \
         'neumann.table1' u 1:4 w lp t 'full cells (J and C, 1998)', \
         'neumann.table1' u 1:2 w lp t 'partial cells (J and C, 1998)'
    Maximum residual convergence (script)

    Maximum residual convergence (script)

    This leads to third-order overall convergence.

    set xrange [*:*]
    fit f(x) '../dirichlet/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:1024]
    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.table2' u 1:2 w lp t 'max all cells (J and C, 1998)',\
         '' u 1:3 w lp t 'avg partial cells (J and C, 1998)', \
         '' u 1:4 w lp t 'avg full cells (J and C, 1998)'
    Error convergence (script)

    Error convergence (script)

    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.

    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:1024]
    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), \
         'neumann.table5' u 1:2 w lp t 'full cells (J and C, 1998)', \
         'neumann.table5' u 1:3 w lp t 'partial cells (J and C, 1998)'
    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($3)) via a,b
    fit f2(x) '' u (log($1)):(log($5)) via a2,b2
    set ylabel 'Maximum error'
    set xrange [16:1024]
    set yrange [*:*]
    plot '' u 1:3 pt 6 t 'all cells', exp(f(log(x))) t ftitle(a,b), \
         '' u 1:5 t 'partial cells', exp(f2(log(x))) t ftitle(a2,b2), \
         'neumann.table5' u 1:5 w lp t 'all cells (J and C, 1998)', \
         'neumann.table5' u 1:4 w lp t 'partial cells (J and C, 1998)'
    Maximum error convergence (script)

    Maximum error convergence (script)

    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 ]