sandbox/Antoonvh/test_poisson_2d.c

    A test for the Poisson solver

    This page concerns the convergence of the Poisson solver for f(x,y),

    \displaystyle \nabla^2 f = s.

    Since sine-rich solutions are not suitable for testing, we design another test case where s is a dipole,

    \displaystyle s = e^{-(x - 2)^2 - y^2} - e^{-(x + 2.)^2 - y^2} 

    #include "grid/multigrid.h" 
#include "poisson.h"

scalar s[], f[];

#define SOURCE (exp(-sq(x - 2.) - sq(y)) - exp(-sq(x + 2.) - sq(y)))

    The cell-averaged source term is approximated with 64 points.

    double source (Point point) {
  double a = 0;
  foreach_child() foreach_child() foreach_child() a += SOURCE;
  return a/pow(pow(2., (double)dimension), 3.);
}

    A reference solution

    Since it may be tricky to find an expression for the cell-averaged values of f, we compute a reference solution using a superior resolution grid. The solution on all relevant levels is then stored in an array.

    int maxlevel = 9;
int reference_level = 11;

double * sol;

void obtain_reference_solution () {
  long c = 0;
  int lr = 0;
  while (lr <= maxlevel)
    c += (1 << (dimension*lr++));
  sol = (double*)malloc(c*sizeof(double));
  init_grid (1 << reference_level);
  foreach() 
    s[] = source (point);
  poisson (f, s);
  restriction ({f});
  long ind = 0;
  for (int l = 0; l <= maxlevel; l++)
    foreach_level(l)
      sol[ind++] = f[];
}

    The setup

    Periodic conditions are set to exclude the error-inducing effects of the approximate values of ghost cells. Furthermore we use a large domain.

    int main() {
  periodic (bottom);
  periodic (left);
  TOLERANCE = 1e-7;
  L0 = 20.;
  X0 = Y0 = -L0/2;
  obtain_reference_solution ();

    We study the convergence starting from a 3-level grid down to the maxlevel.

      for (int l = 3; l <= maxlevel; l++) {
    init_grid (1 << l);
    foreach() 
      s[] = source (point);
    poisson (f, s);

    We skip to the relevant part of the array and compute a total error.

        int lr = 0; long ind = 0;
    while (lr < l)
      ind += (1 << (dimension*lr++));
    double err = 0;
    foreach_level (l) 
      err += fabs(f[] - sol[ind++])*sq(Delta);
    printf ("%d\t%g\n", N, err);
  }
  free (sol);
}

    Result

    We obtain second order-accuracy

    set xr [4 : 1024]
set logscale xy 4
set size square
set grid 
set xlabel 'N'
set ylabel 'Total abs. error'
plot 'out' t 'data', 1000*x**-2 t 'second order'
    (script)

    (script)