hydrostatic2.c

Hydrostatic balance with refined embedded boundaries

Hydrostatic balance with refined embedded boundaries

This test case is related to the test case 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 "grid/quadtree.h"
#include "../myembed.h"
#include "../mycentered.h"
#include "view.h"

Geometry

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

void p_shape (scalar c, face vector f)
{
  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, c, f);
  fractions_cleanup (c, f,
		     smin = 1.e-14, cmin = 1.e-14);
}

Setup

We define the mesh adaptation parameters.

#define lmin (6) // Min mesh refinement level
#define lmax (8) // Max mesh refinement level

int main()
{

The domain is 1\times 1.

  L0 = 1.;
  size (L0);
  origin (-L0/2., -L0/2.);

We set the maximum timestep.

  DT = 1.e-2;

We set the tolerance of the Poisson solver.

  TOLERANCE = 1.e-6;

We initialize the grid.

  N = 1 << (lmin);
  init_grid (N);
  
  run();
}

Boundary conditions

Properties

Initial conditions

event init (i = 0)
{

We define the gravity acceleration vector.

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

We use “third-order” face flux interpolation.

#if ORDER2
  for (scalar s in {u, p, pf})
    s.third = false;
#else
  for (scalar s in {u, p, pf})
    s.third = true;
#endif // ORDER2
  
#if TREE

When using TREE and in the presence of embedded boundaries, we should also define the gradient of u at the cell center of cut-cells.

#endif // TREE

We initialize the embedded boundary.

#if TREE

When using TREE, we refine the mesh around the embedded boundary.

  astats ss;
  int ic = 0;
  do {
    ic++;
    p_shape (cs, fs);
    ss = adapt_wavelet ({cs}, (double[]) {1.e-2},
			maxlevel = (lmax), minlevel = (1));
  } while ((ss.nf || ss.nc) && ic < 100);
#endif // TREE
    
  p_shape (cs, fs);

We also define the volume fraction at the previous timestep csm1=cs.

  csm1 = cs;

We define the no-slip boundary conditions for the velocity.

  u.n[embed] = dirichlet (0);
  u.t[embed] = dirichlet (0);
  p[embed]   = neumann (0);

  uf.n[embed] = dirichlet (0);
  uf.t[embed] = dirichlet (0);
  pf[embed]   = neumann (0);

We finally give an initial guess for the pressure.

  foreach() {
    p[] = cs[] ? (g.x[]*x + g.y[]*y) : nodata; // exact pressure
    pf[] = p[];
  }
  boundary ((scalar *) {p, pf});
}

Embedded boundaries

Outputs

We only solve the Euler equations for a few time iterations.

event logfile (i++; i <= 100) 
{

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

  assert (normf(u.x).max < 1.e-14 &&
	  normf(u.y).max < 1.e-14);
  
  fprintf (stderr, "%d %g %g %d %d %d %d %d %d %g %g %g %g %g %g\n",
	   i, t, dt,
	   mgpf.i, mgpf.nrelax, mgpf.minlevel,
	   mgp.i,  mgp.nrelax,  mgp.minlevel,
	   mgpf.resb, mgpf.resa,
	   mgp.resb,  mgp.resa,
	   normf(u.x).max, normf(u.y).max
	   );
}

Snapshots

event snapshot (t = end)
{
  view (fov = 20,
	tx = 0., ty = 0.,
	bg = {1,1,1},
	width = 800, height = 800);
    
  draw_vof ("cs", filled = -1, lw = 5, fc = {1,1,1});
  squares ("p", spread = -1);
  save ("p.png");
}

Results

The pressure is hydrostatic, in each of the pores.