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 the same porous medium as in porous.c but with 200 random disks so that the pores are entirely closed.

void porous (scalar cs, face vector fs)
{
  int ns = 200; // 160, 80
  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 φ[];
  foreach_vertex() {
    φ[] = HUGE;
    for (double xp = -L0; xp <= L0; xp += L0)
      for (double yp = -L0; yp <= L0; yp += L0)
	for (int i = 0; i < ns; i++)
	  φ[] = intersection (φ[], (sq(x + xp - xc[i]) +
					sq(y + yp - yc[i]) - sq(R[i])));
  }
  fractions (φ, 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 ("defaults");
  event ("metric");

  porous (cs, fs);
  refine (level < 8 && cs[] > 0 && cs[] < 1);
  porous (cs, fs);

  const face vector G[] = {1.,1.};
  a = G;

The system is quite stiff.

  TOLERANCE = 1e-6;
  NITERMAX = 100;
  mgp.nrelax = 100;
  α = 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[] -= α.x[] ? dt*α.x[]*face_gradient_x (p, 0) : 0.;
  boundary ((scalar *){uf});
  
  face vector gf[];
  foreach_face()
    gf.x[] = fm.x[] ? fm.x[]*a.x[] - α.x[]*(p[] - p[-1])/Δ : 0.;
  boundary_flux ({gf});
  
  trash ({g});
  foreach()
    foreach_dimension()
      g.x[] = (gf.x[] + gf.x[1])/(fm.x[] + fm.x[1] + SEPS);
  boundary ((scalar *){g});

  correction (dt);
#endif

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

  fprintf (stderr, "mgp %g %g %d %d\n", mgp.resb, mgp.resa, mgp.i, mgp.minlevel);
  fprintf (stderr, "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");

  scalar p1[];
  foreach()
    p1[] = p[];
  dump();
}