poisson_f.c

Errors for the Poisson solver.
A reference solution
- The setup
- Results

Errors 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 an other test case where s is a dipole,

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

The solution (f) only exists at the mercy of boundary conditions for f.

#include "utils.h"
#include "poisson.h"
scalar s[], f[];

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

f[top] = dirichlet (0.);
f[bottom] = dirichlet (0.);

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

It may be hard to find an analytical expression for the cell-averaged solution f and hence we compute a reference solution via a superior-resolution grid. The solution on all relevant levels is then stored in an array.

#define OFFSET(l) ((1 << (2*(l)))*((l) + 1))
#define _O (-GHOSTS)
#define C_IND(i,j,l) ((i+_O) + (1 << l)*(j+_O))
#define INDEX (OFFSET(level - 1) + C_IND(point.i, point.j, level)) 

int maxlevel = 9;         //The maximum resolution for analysis
int reference_level = 11;  //The superior resolution for reference 

double * sol;

void obtain_reference_solution() {
  sol = (double*)malloc (OFFSET (maxlevel)*sizeof(double));
  init_grid (1 << reference_level);
  foreach() 
    s[] = source (point);
  poisson (f, s);
  restriction ({f});
  foreach_cell()
    if (level <= maxlevel)
      sol[INDEX] = f[];
  output_ppm (s, file = "s.png", n = 350);
  output_ppm (f, file = "f.png", n = 350);
}

The test system looks like this:

The source term (s)

The Solution (f)

The setup

We set a large domain, a small TOLERANCE and a periodic solution in the left-right direction. Then a reference solution is computed.

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

Solutions are approximated on increasingly refined grids. Furthermore, some error data is written.

  FILE * fp = fopen ("error_data", "w");
  for (int l = 4; l <= maxlevel; l++) {
    init_grid (1 << l);
    foreach() 
      s[] = source (point);
    poisson (f, s);
    boundary ({f});
    scalar dfdx[];
    foreach()
      dfdx[] = (f[1] - f[-1])/(2*Delta);
    boundary ({dfdx});
    double err = 0;
    foreach() { 
      double e = f[] - sol[INDEX];
      err += fabs(e)*sq(Delta);
      fprintf (fp, "%d\t%g\t%g\t%g\n", l, e/sq(Delta), s[],
	       (dfdx[0,1] - dfdx[0,-1])/(2*Delta));
    }
    printf ("%d\t%g\n", N, err);
    if (l == maxlevel) {
      output_ppm (s, file = "s.png", n = 512);
      output_ppm (f, file = "f.png", n = 512);
    }
  }
  fclose (fp);
  free (sol);
}

Results

We check the converegence by plotting the total error versus N.

set xr [10:800]
set logscale xy
set xlabel 'N'
set ylabel 'Total error' 
plot 'out' t 'data', 100*x**(-2) t 'second order convergence'

It’s seconder order (script)

This motivates to plot the scaled error against the cell averaged value of s (s[]).

reset
set term pngcairo
set output 'plot1.png'
set xlabel 'e/sq(Delta)'
set ylabel 'level + s[]' 
set grid
plot 'error_data' u 2:($1+$3):1 palette t 'data', -x*16.5 + 7 lw 2 t 'linear'

Error behaviour amongst levels (script)

The error could maybe also scale with \frac{\mathrm{d}^2f}{\mathrm{d}x\mathrm{d}y},

reset
set term pngcairo
set output 'plot2.png'
set xlabel 'e/sq(Delta)'
set ylabel 'level + ddf/dxdy[]' 
set grid
plot 'error_data' u 2:($1+$4):1 palette t 'data'

Error behaviour amongst levels (script)

but is does not.