sandbox/Antoonvh/pb.c

    A dirichlet BC trick for the Poisson solver

    We solve \nabla^2 f = s for f and want to prescribe some value at the top boundary. Normally it would look like this,

    #include "poisson.h"
#include "utils.h"

scalar f[], s[];

int main() {
  periodic (left);
  L0 = 10;
  X0 = -5.;
  Y0 = -2.5;
  init_grid (128);
  foreach()
    s[] = exp(-sq(x) - sq(y));
  output_ppm (s, file = "s.png", n = 256, linear = true);
    The source term looks like this

    The source term looks like this

    We set boundary conditions and solve the problem

      f[top] =  dirichlet (1.);
  f[bottom] = dirichlet (0.);
  poisson (f, s);
  output_ppm (f, file = "f.png", n = 256,
	      linear = true, min = -2, max = 1.5);
    The solution for f on the square domain

    The solution for f on the square domain

    Boundary within the domain

    We want to prescribe the value of f halfway the domain (i.e. y = Y0 + L0/2) and for some reason we do not use mask or embedded boundaries. We must then first set a homogeneous dirichlet condition at the top, so that we solve for a modified field f_{\mathrm{mod}} = f-1.

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

    The dirichlet condition is extented into the domain via a face vector field and a scalar field:

      face vector as[];
  scalar ls[];
# define OUTSIDE_DOMAIN (y >= (Y0 + L0/2.)) 
  foreach() {
    if (OUTSIDE_DOMAIN) { // Outide the internal domain:
      s[] = 0;            // Source term is zero
      ls[] = -10;         // Linear term 
    }else
      ls[] = 0;           // No linear term
  }
  foreach_face()          // Set the default face weights 
    as.x[] = 1;
  foreach_face(x)         // Except for the boundary-perpendicular component,
    if (OUTSIDE_DOMAIN)   // outside the internal domain.
      as.x[] = 0.;
  boundary((scalar*) all);

    We now solve the modified problem:

      poisson (f, s, as, ls);
  output_ppm (f, file = "f2.png", n = 256,
	      linear = true, min = -3, max = 0.5);
    The solution looks as if it could be true

    The solution looks as if it could be true

    we compare it against the proper mask implementation,

      mask (OUTSIDE_DOMAIN ? top : none);
  boundary(all);
  poisson (f, s);
  output_ppm (f, file = "f3.png", n = 256,
	      linear = true, min = -3, max = 0.5);
    Its OK-ish

    Its OK-ish

    }