sandbox/ecipriano/test/poisson-gfm.c

    Test cases for the Poisson–Helmholtz solver with Ghost Fluid Method

    We reproduce the test cases proposed by Liu et al. 2000 to verify the implementation of the GFM approach to impose interface boundary conditions in the solution of Poisson–Helmholtz equations.

    #include "grid/multigrid1D.h"
    #include "poisson-gfm.h"
    #include "fractions.h"
    #include "run.h"
    
    int maxlevel = 11;
    
    scalar a[], f[];
    
    int main (void) {
      L0 = 1 [0];
      init_grid (1 << maxlevel);

    Case 1: Laplace equation

    We solve the following Laplace equation: \displaystyle \nabla^2 a = 0

    with: \displaystyle \left[a\right]_\Gamma = 1 \displaystyle \left[\nabla a\cdot\mathbf{n}_\Gamma\right]_\Gamma = 0

      fraction (f, x - 0.5*L0);
    
      a[left] = dirichlet (0.);
      a[right] = dirichlet (2.);
    
      scalar b[];
      foreach()
        b[] = 0.;
    
      face vector ajump[];
      foreach_face()
          ajump.x[] = 1.;
    
      poisson_ghost (a, b, f = f, ajump = ajump,
          tolerance = TOLERANCE/sq(dt), nrelax = 4);
    
      for (double x = 0.; x < L0; x += 0.5*L0/(1 << maxlevel))
        fprintf (stderr, "%g %g %g\n", x, interpolate (f, x), interpolate (a, x));
      fprintf (stderr, "\n\n");
    reset
    set xr[0:1]
    set grid
    set xlabel "x"
    set ylabel "a"
    set key top left
    
    f(x) = (x <= 0.5) ? x : x + 1
    
    plot "log" index 0 u 1:3 w l lw 2 t "Numerical", \
         f(x) w p pt 34 t "Exact"
    (script)

    (script)

    Case 2: Laplace equation

    We solve the following Laplace equation: \displaystyle \nabla^2 a = 0

    with: \displaystyle \left[a\right]_\Gamma = 0 \displaystyle \left[\nabla a\cdot\mathbf{n}_\Gamma\right]_\Gamma = 1

      a[left] = dirichlet (0.);
      a[right] = dirichlet (1.5);
    
      foreach() {
        a[] = 0.;
        b[] = 0.;
      }
    
      face vector bjump[];
      foreach_face()
        bjump.x[] = 1.;
    
      poisson_ghost (a, b, f = f, bjump = bjump,
          tolerance = TOLERANCE/sq(dt), nrelax = 4);
    
      for (double x = 0.; x < L0; x += 0.5*L0/(1 << maxlevel))
        fprintf (stderr, "%g %g %g\n", x, interpolate (f, x), interpolate (a, x));
      fprintf (stderr, "\n\n");
    reset
    set xr[0:1]
    set grid
    set xlabel "x"
    set ylabel "a"
    set key top left
    
    f(x) = (x <= 0.5) ? x : 2.*x - 0.5
    
    plot "log" index 1 u 1:3 w l lw 2 t "Numerical", \
         f(x) w p pt 34 t "Exact"
    (script)

    (script)

    Case 3: Poisson equation

    We solve the following Poisson equation: \displaystyle \nabla\cdot\left(\beta\nabla a\right) = \begin{cases} -2e^x\sin (x) & \text{ if } x < 1 \\ 4(\cos^2(x) - \sin^2(x)) & \text { if } x \geq 1 \end{cases}

    with: \displaystyle \left[a\right]_\Gamma = e^x\cos (x) - (\sin^2(x) - \cos^2(x)) \displaystyle \left[\beta\nabla a\cdot\mathbf{n}_\Gamma\right]_\Gamma = - e^x(\cos (x) - \sin (x)) - 4\sin (x) \cos (x)

      L0 = 2 [0];
      init_grid (1 << maxlevel);
    
      a[left] = dirichlet (exp(x)*cos(x));
      a[right] = dirichlet (sq(sin(x)) - sq(cos(x)));
    
      foreach()
        f[] = (x < 1.) ? 1. : 0.;
    
      foreach() {
        a[] = 0.;
        b[] = (x < 1.) ? (-2.*exp(x)*sin(x)) : (4.*(sq(cos(x)) - sq(sin(x))));
      }
    
      face vector betac[];
      foreach_face()
        betac.x[] = 1.;
    
      foreach_face()
        ajump.x[] = (exp(x)*cos(x) - sq(sin(x)) + sq(cos(x)));
    
      foreach_face()
        bjump.x[] = -(exp(x)*(cos(x) - sin(x)) - 4.*sin(x)*cos(x));
    
      poisson_ghost (a, b, f = f, ajump = ajump, bjump = bjump, alpha = betac,
          tolerance = TOLERANCE/sq(dt), nrelax = 4);
    
      for (double x = 0.; x < L0; x += 0.5*L0/(1 << maxlevel))
        fprintf (stderr, "%g %g %g\n", x, interpolate (f, x), interpolate (a, x));
      fprintf (stderr, "\n\n");
    reset
    set xr[0:2]
    set grid
    set xlabel "x"
    set ylabel "a"
    set key top left
    
    f(x) = (x <= 1.) ? (exp(x)*cos(x)) : (sin(x)**2 - cos(x)**2)
    
    plot "log" index 2 u 1:3 w l lw 2 t "Numerical", \
         f(x) w p pt 34 t "Exact"
    (script)

    (script)

    }

    References

    [liu2017second]

    Cheng Liu and Changhong Hu. A second order ghost fluid method for an interface problem of the poisson equation. Communications in Computational Physics, 22(4):965–996, 2017.

    [liu2000boundary]

    Xu-Dong Liu, Ronald P Fedkiw, and Myungjoo Kang. A boundary condition capturing method for poisson’s equation on irregular domains. Journal of computational Physics, 160(1):151–178, 2000.