src/test/hydrostatic2.c

    Hydrostatic balance with refined embedded boundaries

    This test case is related to flow in a complex porous medium. We check that for a “closed pore” medium, hydrostatic balance can be recovered. This is not trivial in particular when the spatial resolution is variable, since pressure values need to be interpolated (at least) to second-order close to embedded boundaries. This behaviour depends on the restriction and prolongation operators used close to embedded boundaries.

    #include "embed.h"
    #include "navier-stokes/centered.h"
    #include "view.h"

    We use a similar porous medium as in porous.c but with enough disks so that the pores are entirely closed.

    void porous (scalar cs, face vector fs)
    {
      int ns = 200;
      double xc[ns], yc[ns], R[ns];
      srand (0);
      for (int i = 0; i < ns; i++)
        xc[i] = 0.5*noise(), yc[i] = 0.5*noise(), R[i] = 0.02 + 0.04*fabs(noise());
    
      vertex scalar phi[];
      foreach_vertex() {
        phi[] = HUGE;
        for (double xp = -L0; xp <= L0; xp += L0)
          for (double yp = -L0; yp <= L0; yp += L0)
    	for (int i = 0; i < ns; i++)
    	  phi[] = intersection (phi[], (sq(x + xp - xc[i]) +
    					sq(y + yp - yc[i]) - sq(R[i])));
      }
      boundary ({phi});
      fractions (phi, cs, fs);
      fractions_cleanup (cs, fs);
    }
    
    int main()
    {
      origin (-0.5, -0.5);
    
      init_grid (1 << 6);

    The events of the Navier-Stokes solver are called “by hand”.

      event ("metric");
      event ("defaults");
    
      porous (cs, fs);
    #if 0
      refine (level < 8 && cs[] > 0 && cs[] < 1);
      porous (cs, fs);
    #else
      adapt_wavelet ({cs}, (double[]){0.01}, 8, 6);
      porous (cs, fs);
      adapt_wavelet ({cs}, (double[]){0.01}, 8, 6);
      porous (cs, fs);
      adapt_wavelet ({cs}, (double[]){0.01}, 8, 6);
      porous (cs, fs);
    #endif
      
      const face vector G[] = {1.,2.};
      a = G;

    The system is quite stiff.

      TOLERANCE = 1e-6;
      NITERMAX = 100;
      mgp.nrelax = 100;
      alpha = fm;
      dt = 1.;
    
      event ("acceleration");
    #if 1
      event ("projection");
    #else
      foreach()
        p[] = G.x[]*x + G.y[]*y; // exact pressure
      boundary ({p});
      foreach_face()
        uf.x[] -= alpha.x[] ? dt*alpha.x[]*face_gradient_x (p, 0) : 0.;
      boundary ((scalar *){uf});
      
      face vector gf[];
      foreach_face()
        gf.x[] = fm.x[] ? fm.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta : 0.;
      boundary_flux ({gf});
      
      trash ({g});
      foreach()
        foreach_dimension() {
          g.x[] = (gf.x[] + gf.x[1])/(fm.x[] + fm.x[1] + SEPS);
          //      assert (fabs(g.x[]) < 1e-6);
        }
      boundary ((scalar *){g});
    
      correction (dt);
    #endif

    We check the convergence rate and the norms of the velocity field (which should be negligible).

      fprintf (ferr, "mgp %g %g %d %d\n", mgp.resb, mgp.resa, mgp.i, mgp.minlevel);
      fprintf (ferr, "umax %g %g\n", normf(u.x).max, normf(u.y).max);

    The pressure is hydrostatic, in each of the pores.

    Pressure field.

    Pressure field.

      view (fov = 19, width = 400, height = 400);
      draw_vof ("cs", filled = -1, fc = {1,1,1});
      squares ("p", spread = -1);
      save ("p.png");
    
      p.nodump = false;
      scalar pid[];
      foreach()
        pid[] = pid();
      dump();
    }