src/embed.h

Embedded boundaries

Boundaries of general shape can be described using an integral (i.e. finite volume) formulation which takes into account the volume and area fractions of intersection of the embedded boundary with the Cartesian mesh.

We will need to deal with volume fractions. Interpolations (for Dirichlet boundary conditions) assume a 5x5 stencil.

#include "fractions.h"
#define BGHOSTS 2
#define EMBED 1

The volume and area fractions are stored in these fields.

scalar cs[];
face vector fs[];

Operator overloading

Several standard operators, defined in common.h need to be tuned to take into account the embedded fractions.

The SEPS constant is used to avoid division by zero.

#undef SEPS
#define SEPS 1e-30

Face gradients and face values, computed from cell-centered values must be tuned to take into account the area fractions of the embedded boundary. This follows the procedure described in Johansen and Colella, 1998, figure 3 and equation (16) in particular.

foreach_dimension()
static inline double embed_face_gradient_x (Point point, scalar a, int i)
{
  assert (dimension == 2);
  int j = sign(fs.x[i,1] - fs.x[i,-1]);
  return fs.x[i,j] > 0.5 ?
    ((1. + fs.x[i])*(a[i] - a[i-1]) +
     (1. - fs.x[i])*(a[i,j] - a[i-1,j]))/(2.*Δ) :
    (a[i] - a[i-1])/Δ;
}

foreach_dimension()
static inline double embed_face_value_x (Point point, scalar a, int i)
{
  assert (dimension == 2);
  int j = sign(fs.x[i,1] - fs.x[i,-1]);
  return fs.x[i,j] > 0.5 ?
    ((1. + fs.x[i])*(a[i]*cs[i] + a[i-1]*cs[i-1])/(cs[i] + cs[i-1] + SEPS) +
     (1. - fs.x[i])*(a[i,j]*cs[i,j] + a[i-1,j]*cs[i-1,j])/
     (cs[i,j] + cs[i-1,j] + SEPS))/2. :
    (a[i]*cs[i] + a[i-1]*cs[i-1])/(cs[i] + cs[i-1] + SEPS);
}

We use the functions above to redefine the face gradient macros.

#undef face_gradient_x
#define face_gradient_x(a, i)						\
  (fs.x[i] > 0. && fs.x[i] < 1. ?					\
   embed_face_gradient_x (point, a, i) :				\
   (a[i] - a[i-1])/Δ)

#undef face_gradient_y
#define face_gradient_y(a, i)						\
  (fs.y[0,i] > 0. && fs.y[0,i] < 1. ?					\
   embed_face_gradient_y (point, a, i) :				\
   (a[0,i] - a[0,i-1])/Δ)

#undef face_gradient_z
#define face_gradient_z(a, i)						\
  (fs.z[0,0,i] > 0. && fs.z[0,0,i] < 1. ?				\
   embed_face_gradient_z (point, a, i) :				\
   (a[0,0,i] - a[0,0,i-1])/Δ)

The face values, interpolated from cell-centered values, need to use weighting by the cell fraction for stability (of the centered Navier–Stokes solver). The corresponding test case is test/uf.c.

#undef face_value
#define face_value(a, i)					\
  (fs.x[i] > 0. && fs.x[i] < 1. ?				\
   embed_face_value_x (point, a, i) :				\
   (a[i]*cs[i] + a[i-1]*cs[i-1])/(cs[i] + cs[i-1] + SEPS))

The centered gradient must not use values of fields entirely contained within the embedded boundary (for which cs is zero).

#undef center_gradient
#define center_gradient(a) (cs[1] && cs[-1] ? (a[1] - a[-1])/(2.*Δ) : \
			    cs[1] && cs[]   ? (a[1] - a[])/Δ :	  \
			    cs[-1] && cs[]  ? (a[] - a[-1])/Δ : 0.)

Utility functions for the geometry of embedded boundaries

For a cell containing a fragment of embedded boundary (i.e. for which 0<cs<1), embed_geometry() returns the area of the fragment, the relative position p of the barycenter of the fragment and the boundary normal n.

static inline
double embed_geometry (Point point, coord * p, coord * n)
{
  *n = facet_normal (point, cs, fs);
  double α = plane_alpha (cs[], *n);
  double area = plane_area_center (*n, α, p);
  normalize (n);
  return area;
}

This function and the macro below shift the position (x1,y1,z1) to the position of the barycenter of the embedded fragment.

static inline
double embed_area_center (Point point, double * x1, double * y1, double * z1)
{
  double area = 0.;
  if (cs[] > 0. && cs[] < 1.) {
    coord n, p;
    area = embed_geometry (point, &p, &n);
    *x1 += p.x*Δ, *y1 += p.y*Δ, *z1 += p.z*Δ;
  }
  return area;
}

#define embed_pos() embed_area_center (point, &x, &y, &z)

This function returns the value of field s interpolated linearly at the barycenter p of the fragment of embedded boundary contained within the cell.

double embed_interpolate (Point point, scalar s, coord p)
{
  assert (dimension == 2);
  int i = sign(p.x), j = sign(p.y);
  if (cs[i] && cs[0,j] && cs[i,j])
    // bilinear interpolation when all neighbors are defined
    return ((s[]*(1. - fabs(p.x)) + s[i]*fabs(p.x))*(1. - fabs(p.y)) + 
	    (s[0,j]*(1. - fabs(p.x)) + s[i,j]*fabs(p.x))*fabs(p.y));
  else {
    // linear interpolation with gradients biased toward the
    // cells which are defined
    double val = s[];
    foreach_dimension() {
      int i = sign(p.x);
      if (cs[i])
	val += fabs(p.x)*(s[i] - s[]);
      else if (cs[-i])
	val += fabs(p.x)*(s[] - s[-i]);
    }
    return val;
  }
}

Dirichlet boundary condition

This function returns the gradient of scalar s, normal to the embedded boundary defined by cs, of unit normal vector n (normalised using the Euclidean norm, not the box norm) and of centroid p. The Dirichlet boundary condition bc is imposed on the embedded boundary.

The calculation follows Johansen and Colella, 1998 and is summarised in the figure below (see also Figure 4 of Johansen and Colella).

Third-order normal gradient scheme

Third-order normal gradient scheme

Note that the function will return nodata if the gradient cannot be computed (which can happen for degenerate geometric configurations).

double dirichlet_gradient (Point point, scalar s, scalar cs,
			   coord n, coord p, double bc)
{
  assert (dimension == 2);
  n.x = - n.x, n.y = - n.y;
  foreach_dimension()
    if (fabs(n.x) >= fabs(n.y)) {
      double d[2], v[2];
      for (int k = 0; k <= 1; k++) {
	int i = (k + 1)*sign(n.x);
	d[k] = (i - p.x)/n.x;
	double y = p.y + d[k]*n.y;
	int j = y > 0.5 ? 1 : y < -0.5 ? -1 : 0;
	y -= j;
	if (cs[i,j] && cs[i,j - sign(n.y)] && cs[i,j + sign(n.y)] > 0.5)
	  v[k] = (s[i,j-1]*(y - 1.) + s[i,j+1]*(y + 1.))*y/2.
	    - s[i,j]*(y + 1.)*(y - 1.); // third-order interpolation
	else
	  v[k] = nodata;
      }
      if (v[0] != nodata) {
	if (v[1] != nodata) // third-order gradient
	  return (d[1]*(bc - v[0])/d[0] - d[0]*(bc - v[1])/d[1])/
	    ((d[1] - d[0])*Δ);
	return (bc - v[0])/(d[0]*Δ); // second-order gradient
      }
    }
  return nodata; // degenerate case
}

bid embed;

Surface force and vorticity

We first define a function which computes \nablaun while taking the boundary conditions on the embedded surface into account.

static inline
coord embed_gradient (Point point, vector u, coord p, coord n)
{
  coord dudn;
  foreach_dimension() {
    double sb = u.x[], vb = u.x.boundary[embed] (point, point, u.x);
    bool dirichlet = u.x[]; u.x[] = sb;
    if (dirichlet)
      dudn.x = dirichlet_gradient (point, u.x, cs, n, p, vb);
    else // Neumann
      dudn.x = vb;
    if (dudn.x == nodata)
      dudn.x = 0.;
  }
  return dudn;
}

The force exerted by the fluid on the solid can be written FΓ=Γ(pI+2μD)ndΓ with Γ the solid boundary. It can be further decomposed into a pressure (i.e. “form”) drag Fp=ΓpndΓ and a viscous drag Fμ=Γ2μDndΓ These two vectors are computed by the embed_force() function.

trace
void embed_force (scalar p, vector u, face vector mu, coord * Fp, coord * Fmu)
{
  // fixme: this could be simplified considerably if reduction worked on vectors
  double Fpx = 0., Fpy = 0., Fmux = 0., Fmuy = 0.;
  foreach (reduction(+:Fpx) reduction(+:Fpy)
	   reduction(+:Fmux) reduction(+:Fmuy))
    if (cs[] > 0. && cs[] < 1.) {

To compute the pressure force, we first get the coordinates of the barycentre of the embedded fragment, its area and normal, and then interpolate the pressure field on the surface.

      coord n, b;
      double area = embed_geometry (point, &b, &n);
      area *= pow (Δ, dimension - 1);
      double Fn = area*embed_interpolate (point, p, b);
      Fpx += Fn*n.x;
      Fpy += Fn*n.y;

To compute the viscous force, we first need to retrieve the local value of the viscosity (ideally at the barycentre of the embedded fragment). This is not completely trivial since it is defined on the faces of the cell. We use a surface-fraction-weighted average value.

      if (constant(μ.x) != 0.) {
	double mua = 0., fa = 0.;
	foreach_dimension() {
	  mua += μ.x[] + μ.x[1];
	  fa  += fs.x[] + fs.x[1];
	}
	mua /= fa;

To compute the viscous force, we need to take into account the (Dirichlet) boundary conditions for the velocity on the surface. We only know how to do this when computing the normal gradient \nablaun using the embed_gradient() function. We thus need to re-express the viscous force using only normal derivatives of the velocity field.

If we assume that u is constant on the boundary, then $ \mathbf{{\nabla}} \mathbf{u} \cdot \mathbf{t}= \mathbf{0} $ with t the unit tangent vector to the boundary. We thus have the relations $ \mathbf{{\nabla}} \mathbf{u} = \left( \mathbf{{\nabla}} \mathbf{u} \cdot \mathbf{n} \right) \mathbf{n} + \left( \mathbf{{\nabla}} \mathbf{u} \cdot \mathbf{t} \right) \mathbf{t} = \left( \mathbf{{\nabla}} \mathbf{u} \cdot \mathbf{n} \right) \mathbf{n} $ $ \mathbf{D}= \frac{1}{2} \left( \mathbf{{\nabla}} \mathbf{u} + \mathbf{{\nabla}}^T \mathbf{u} \right) = \frac{1}{2} \left(\begin{array}{cc} 2 \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_x & \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x\\ \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x & 2 \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_y \end{array}\right) $ $ \mathbf{F}_{\mu} = - \int_{\Gamma} \left(\begin{array}{c} \left[2 \mu \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_x \right] n_x + \mu \left[ \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x \right] n_y\\ \left[2 \mu \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_y \right] n_y + \mu \left[ \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x \right] n_x \end{array}\right) $ $ \mathbf{F}_{\mu} = - \int_{\Gamma} \left(\begin{array}{c} \mu \left[ \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) (n^2_x + 1) + \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x n_y \right]\\ \mu \left[ \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) (n^2_y + 1) + \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_x n_y \right] \end{array}\right) $

	assert (dimension == 2);
	coord dudn = embed_gradient (point, u, b, n);
	Fmux -= area*mua*(dudn.x*(sq(n.x) + 1.) + dudn.y*n.x*n.y);
	Fmuy -= area*mua*(dudn.y*(sq(n.y) + 1.) + dudn.x*n.x*n.y);
      }
    }

  Fp->x = Fpx; Fp->y = Fpy; 
  Fmu->x = Fmux; Fmu->y = Fmuy; 
}

In 2D dimensions, embed_vorticity() returns the vorticity of velocity field u, on the surface of the embedded boundary contained in the cell. p is the relative position of the barycentre of the embedded fragment and n its normal.

#if dimension == 2
double embed_vorticity (Point point, vector u, coord p, coord n)
{

We compute $\mathbf{{\nabla}}\mathbf{u}\cdot\mathbf{n}$, taking the boundary conditions into account.

  coord dudn = embed_gradient (point, u, p, n);

The vorticity is then obtained using the relations $ \omega = \partial_x v - \partial_y u = \left( \mathbf{{\nabla}} v \cdot \mathbf{n} \right) n_x - \left( \mathbf{{\nabla}} u \cdot \mathbf{n} \right) n_y $

  return dudn.y*n.x - dudn.x*n.y;
}
#endif // dimension == 2

Flux through the embedded boundary

This function computes the flux through the embedded boundary contained within a cell bμsndb with db the elementary boundary surface and n the embedded boundary (outward-pointing) normal.

Boundary conditions for s are taken into account (and are considered homogeneous if set to true).

The result is returned in val.

The function returns false if the flux cannot be computed, in which case val is set to the value of the Dirichlet boundary condition.

bool embed_flux (Point point, scalar s, face vector mu,
		 bool homogeneous, double * val)
{

If the cell does not contain a fragment of embedded boundary, the flux is zero.

  *val = 0.;
  if (cs[] >= 1. || cs[] <= 0.)
    return true;

We compute the normal and the barycenter of the fragment of embedded boundary contained within the cell.

  coord n = facet_normal (point, cs, fs), p;
  double α = plane_alpha (cs[], n);
  double area = plane_area_center (n, α, &p);

  double sb = s[], vb = s.boundary[embed] (point, point, s);
  bool dirichlet = s[]; s[] = sb;

If the boundary condition is homogeneous Neumann, the flux is zero.

  if (homogeneous && !dirichlet)
    return true;

If the boundary condition is Dirichlet, we need to compute the normal gradient.

  if (dirichlet) {
    if (homogeneous)
      vb = 0.;
    normalize (&n);
    double dg = dirichlet_gradient (point, s, cs, n, p, vb);

If the gradient cannot be computed, we return the value of the Dirichlet condition.

    if (dg == nodata) {
      *val = vb;
      return false;
    }
    vb = dg;
  }

We retrieve the (average) value of μ without the metric.

  double mua = 0., fa = 0.;
  foreach_dimension() {
    mua += μ.x[] + μ.x[1];
    fa  += fs.x[] + fs.x[1];
  }
  *val = - mua/(fa + SEPS)*vb*area/Δ;
  return true;
}

#define neumann_embed(expr)   (embed_area_center (point, &x, &y, &z), \
			       val(_s) = 0, (expr))
#define dirichlet_embed(expr) (embed_area_center (point, &x, &y, &z), \
			       val(_s) = 1, (expr))

Restriction and refinement/prolongation

Volume fraction field cs

For the embedded fraction field cs, the function below is modelled closely on the volume fraction refinement function fraction_refine().

#if TREE
static void embed_fraction_refine (Point point, scalar cs)
{
  double cc = cs[];

If the cell is empty or full, simple injection from the coarse cell value is used.

  if (cc <= 0. || cc >= 1.) {
    foreach_child()
      cs[] = cc;
  }
  else {

If the cell contains the embedded boundary, we reconstruct the boundary using VOF linear reconstruction and a normal estimated from the surface fractions.

    coord n = facet_normal (point, cs, fs);
    double α = plane_alpha (cc, n);
      
    foreach_child() {
      static const coord a = {0.,0.,0.}, b = {.5,.5,.5};
      coord nc;
      foreach_dimension()
	nc.x = child.x*n.x;
      cs[] = rectangle_fraction (nc, α, a, b);
    }
  }
}

Surface fractions field fs

The embedded surface fractions fs are reconstructed using this function.

foreach_dimension()
static void embed_face_fraction_refine_x (Point point, scalar s)
{
  vector fs = s.v;

If the cell is empty or full, simple injection from the coarse cell value is used.

  if (cs[] <= 0. || cs[] >= 1.) {

We need to make sure that the fine cells face fractions match those of their neighbours.

    fine(fs.x,1,0) = fine(fs.x,1,1) = cs[];
    if (!is_refined(neighbor(-1)) &&
	(is_local(cell) || is_local(neighbor(-1))))
      fine(fs.x,0,0) = fine(fs.x,0,1) = fs.x[];
    if (!is_refined(neighbor(1)) && neighbor(1).neighbors &&
	(is_local(cell) || is_local(neighbor(1))))
      fine(fs.x,2,0) = fine(fs.x,2,1) = fs.x[1];
  }
  else {

If the cell contains the embedded boundary, we reconstruct the boundary using VOF linear reconstruction and a normal estimated from the surface fractions.

    coord n = facet_normal (point, cs, fs);
    double α = plane_alpha (cs[], n);

We need to reconstruct the face fractions fs for the fine cells.

#if dimension == 2

For the fine face fractions contained within the coarse cell, we compute the intersections directly using the VOF reconstruction, and obtain the face fractions by taking into account the orientation of the normal.

    if (2.*fabs(α) < fabs(n.y)) {
      double yc = α/n.y;
      int i = yc > 0.;
      fine(fs.x,1,1 - i) = n.y < 0. ? 1. - i : i;
      fine(fs.x,1,i) = n.y < 0. ? i - 2.*yc : 1. - i + 2.*yc;
    }
    else
      fine(fs.x,1,0) = fine(fs.x,1,1) = α > 0.;

For the fine face fractions coincident with the faces of the coarse cell, we obtain the intersection position from the coarse cell face fraction. The orientation is obtained by looking at the values of face fractions in the transverse direction.

    for (int i = 0; i <= 1; i++)
      if (!is_refined(neighbor(2*i-1)) && neighbor(2*i-1).neighbors &&
	  (is_local(cell) || is_local(neighbor(2*i-1)))) {
	if (fs.x[i] <= 0. || fs.x[i] >= 1.)
	  fine(fs.x,2*i,0) = fine(fs.x,2*i,1) = fs.x[i];
	else {
	  double a = fs.y[0,1] <= 0. || fs.y[2*i-1,1] <= 0. ||
	    fs.y[] >= 1. || fs.y[2*i-1] >= 1.;
	  if ((2.*a - 1)*(fs.x[i] - 0.5) > 0.) {
	    fine(fs.x,2*i,0) = a;
	    fine(fs.x,2*i,1) = 2.*fs.x[i] - a;
	  }
	  else {
	    fine(fs.x,2*i,0) = 2.*fs.x[i] + a - 1.;
	    fine(fs.x,2*i,1) = 1. - a;
	  }
	}
      }
#else // dimension != 2
      assert (false); // not implemented yet
#endif
  }
}

Restriction of cell-centered fields

We now define restriction and prolongation functions for cell-centered fields. The goal is to define second-order operators which do not use any values from cells entirely contained within the embedded boundary (for which cs = 0).

static inline void restriction_embed_linear (Point point, scalar s)
{

For restriction, the cases where 0, 3 or 4 children cells are not entirely contained within the embedded boundary is simple, since linear interpolation can be used directly.

  // 0 children
  if (!cs[]) {
    s[] = 0.;
    return;
  }
  // 4 children
  if (fine(cs,0,0) && fine(cs,1,0) && fine(cs,0,1) && fine(cs,1,1)) {
    s[] = (fine(s,0,0) + fine(s,1,0) + fine(s,0,1) + fine(s,1,1))/4.;
    return;
  }
  // 3 children
  if (fine(cs,0,0) && fine(cs,1,1)) {
    s[] = (fine(s,0,0) + fine(s,1,1))/2.;
    return;
  }
  if (fine(cs,0,1) && fine(cs,1,0)) {
    s[] = (fine(s,0,1) + fine(s,1,0))/2.;
    return;
  }

We are now left with the cases of two children sharing the same coarse face, or with a single child.

In the case of two children, we extrapolate linearly using only the neighboring coarse leaf cell value or fine cell values, as illustrated on the figure below (where the target value is the black dot).

Extrapolation in the case of two children.

Extrapolation in the case of two children.

  // 2 children
  foreach_dimension()
    for (int i = 0; i <= 1; i++)
      if (fine(cs,i,0) && fine(cs,i,1)) {
	if (is_leaf(neighbor(2*i-1)))
	  s[] = cs[2*i-1] ? (2.*(fine(s,i,0) + fine(s,i,1)) - s[2*i-1])/3. :
	    (fine(s,i,0) + fine(s,i,1))/2.;
	else {
	  double v = 0.; // v/n is the average of the fine cell neighbors
	  int n = 0;
	  for (int j = 0; j <= 1; j++)
	    if (fine(cs,3*i-1,j))
	      v += fine(s,3*i-1,j), n++;
	  if (n == 0) {
	    s[] = (fine(s,i,0) + fine(s,i,1))/2.;
	    return;
	  }
	  s[] = (3.*(fine(s,i,0) + fine(s,i,1)) - 2.*v/n)/4.;
	}
	return;
      }

There is a single child not entirely contained within the embedded boundary.

  // 1 child
  for (int i = 0; i <= 1; i++)
    for (int j = 0; j <= 1; j++)
      if (fine(cs,i,j)) {
	if (is_leaf(neighbor(2*i-1,2*j-1))) {
	  assert (cs[2*i-1,2*j-1]);
	  s[] = (4.*fine(s,i,j) - s[2*i-1,2*j-1])/3.;
	}
	else if (fine(cs,3*i-1,3*j-1))
	  s[] = (3.*fine(s,i,j) - fine(s,3*i-1,3*j-1))/2.;
	else
	  s[] = fine(s,i,j);
	return;
      }
  
  assert (false); // not reached
}

Refinement/prolongation of cell-centered fields

For refinement, we use either bilinear interpolation, if the required four coarse cell values are defined or trilinear interpolation if only three coarse cell values are defined. If less that three coarse cell values are defined (“pathological cases” below), we try to estimate gradients in each direction and add the corresponding correction.

static inline void refine_embed_linear (Point point, scalar s)
{
  if (!cs[]) {
    foreach_child()
      s[] = 0.;
    return;
  }
  foreach_child() {
    if (!cs[])
      s[] = 0.;
    else {
      if (coarse(fs.x,(child.x + 1)/2) && coarse(fs.y,0,(child.y + 1)/2) &&
	  coarse(fs.x,(child.x + 1)/2,child.y) &&
	  coarse(fs.y,child.x,(child.y + 1)/2))
	// bilinear interpolation
	s[] = (9.*coarse(s) + 
	       3.*(coarse(s,child.x) + coarse(s,0,child.y)) + 
	       coarse(s,child.x,child.y))/16.;
      else if (coarse(fs.x,(child.x + 1)/2) && coarse(fs.y,0,(child.y + 1)/2))
	// triangular interpolation
	s[] = (2.*coarse(s) + coarse(s,child.x) + coarse(s,0,child.y))/4.;
      else {
	// pathological cases
	s[] = coarse(s);
	foreach_dimension() {
	  if (coarse(fs.x,(child.x + 1)/2))
	    s[] += (coarse(s,child.x) - coarse(s))/4.;
	  else if (!is_leaf(aparent(0,child.y)) &&
		   aparent(0,child.y).neighbors &&
		   cs[0,child.y] && cs[-child.x,child.y])
	    s[] += (s[0,child.y] - s[-child.x,child.y])/2.;
	}
      }
    }
  }
}

Refinement/prolongation of face-centered fields

This function is modelled on refine_face_x() and is typically used to refine the values of the face-centered velocity field uf. It uses linear interpolation, taking into account the weighting by the embedded fractions fs.

void refine_embed_face (Point point, scalar s)
{
  vector v = s.v;
  foreach_dimension() {
    for (int i = 0; i <= 1; i++)
      if (!is_refined(neighbor(2*i - 1)) && neighbor(2*i - 1).neighbors &&
	  (is_local(cell) || is_local(neighbor(2*i - 1)))) {
	double g1 = fs.x[i,+1] && fs.x[i,-1] ?
	  (v.x[i,+1]/fs.x[i,+1] - v.x[i,-1]/fs.x[i,-1])/8. : 0.;
	double g2 = fs.x[i,0,+1] && fs.x[i,0,-1] ?
	  (v.x[i,0,+1]/fs.x[i,0,+1] - v.x[i,0,-1]/fs.x[i,0,-1])/8. : 0.;
	for (int j = 0; j <= 1; j++)
	  for (int k = 0; k <= 1; k++)
	    fine(v.x,2*i,j,k) = fs.x[i] ?
	      fine(fs.x,2*i,j,k)*(v.x[i]/fs.x[i] +
				  (2*j - 1)*g1 + (2*k - 1)*g2) : 0.;
      }
    if (is_local(cell)) {
      double g1 = fs.x[0,+1] && fs.x[1,+1] && fs.x[0,-1] && fs.x[1,-1] ?
	(v.x[0,+1]/fs.x[0,+1] + v.x[1,+1]/fs.x[1,+1] -
	 v.x[0,-1]/fs.x[0,-1] - v.x[1,-1]/fs.x[1,-1])/16. : 0.;
      double g2 = fs.x[0,0,+1] && fs.x[1,0,+1] && fs.x[0,0,-1] && fs.x[1,0,-1] ?
	(v.x[0,0,+1]/fs.x[0,0,+1] + v.x[1,0,+1]/fs.x[1,0,+1] -
	 v.x[0,0,-1]/fs.x[0,0,-1] - v.x[1,0,-1]/fs.x[1,0,-1])/16. : 0.;
      for (int j = 0; j <= 1; j++)
	for (int k = 0; k <= 1; k++)
	  fine(v.x,1,j,k) = fs.x[] && fs.x[1] ?
	    fine(fs.x,1,j,k)*((v.x[]/fs.x[] + v.x[1]/fs.x[1])/2. +
			      (2*j - 1)*g1 + (2*k - 1)*g2) : 0.;
    }
  }    
}

#define refine_face_solenoidal refine_embed_face
#endif // TREE

Prolongation for the multigrid solver

We use a simplified prolongation operator for the multigrid solver i.e. simple injection if bilinear interpolation would use values which are fully contained within the embedded boundary.

#if MULTIGRID
static inline double bilinear_embed (Point point, scalar s)
{
  if (!coarse(cs,0) || !coarse(cs,child.x))
    return coarse(s);
  #if dimension >= 2
  if (!coarse(cs,0,child.y) || !coarse(cs,child.x,child.y))
    return coarse(s);
  #endif
  #if dimension >= 3
  if (!coarse(cs,0,0,child.z) || !coarse(cs,child.x,0,child.z) ||
      !coarse(cs,0,child.y,child.z) ||
      !coarse(cs,child.x,child.y,child.z))
    return coarse(s);  
  #endif
  return bilinear (point, s);
}

#define bilinear(point, s) bilinear_embed(point, s)
#endif // MULTIGRID

Lifting the “small cell” CFL restriction

For explicit advection schemes, the timestep is limited by the CFL conditions Δt<csΔfiui where i is the index of each face, and cs and fi are the embedded volume and face fractions respectively. It is clear that the timestep may need to be arbitrarily small if cs/fs tends toward zero. This is the “small cell” restriction of cut-cell finite-volume techniques.

A classical technique to avoid this limitation is to use a “cell merging” procedure, where the fluxes from cells which would “overflow” are redistributed to neighboring cells.

The function below uses this approach to update a field f, advected by the face velocity field uf, with corresponding advection fluxes flux, during timestep dt which only verifies the standard CFL condition Δt<Δui

trace
void update_tracer (scalar f, face vector uf, face vector flux, double dt)
{

Note that the distinction should be made between cm, the cell fraction metric, and cs, the embedded fraction. This is not done now so that embedded boundaries cannot be combined with a metric yet.

The field e will store the “overflowing” sum of fluxes for each cell.

  scalar e[];
  foreach() {

If the cell is empty, it cannot overflow.

    if (cs[] <= 0.)
      e[] = 0.;

If the cell does not contain an embedded boundary, it cannot overflow either and the sum of the fluxes can be added to advance f in time.

    else if (cs[] >= 1.) {
      foreach_dimension()
	f[] += dt*(flux.x[] - flux.x[1])/Δ;
      e[] = 0.;
    }

If the cell contains the embedded boundary, we compute the maximum timestep verifying the restrictive CFL condition Δtmax=csΔmax(fiui) Note that fs does not appear in the code below because uf already stores the product fsu.

    else {
      double umax = 0.;
      for (int i = 0; i <= 1; i++)
	foreach_dimension()
	  if (fabs(uf.x[i]) > umax)
	    umax = fabs(uf.x[i]);
      double dtmax = Δ*cs[]/(umax + SEPS);

We compute the sum of the fluxes.

      double F = 0.;
      foreach_dimension()
	F += flux.x[] - flux.x[1];
      F /= Δ*cs[];

If the timestep is smaller than Δtmax, the cell cannot overflow and f is advanced in time using the entire flux.

      if (dt <= dtmax) {
	f[] += dt*F;
	e[] = 0.;
      }

Otherwise, the cell is filled “to the brim” by advancing f using the maximum allowable timestep. The field e is used to store the excess flux, weighted by the sum of the neighboring embedded fractions.

      else {
	f[] += dtmax*F;
	double scs = 0.;
	foreach_neighbor(1)
	  scs += sq(cs[]);
	e[] = (dt - dtmax)*F*cs[]/scs;
      }
    }
  }
  boundary ({e});

In a second phase, the excesses in each cell are added to the neighboring cells in proportion of their embedded fractions.

  foreach() {
    double se = 0.;
    foreach_neighbor(1)
      se += e[];
    f[] += cs[]*se;
  }
}

Default settings

To apply the volume/area fraction-weighting to the solvers, we define the metric using the embedded fractions.

event metric (i = 0)
{
  foreach()
    cs[] = 1.;
  foreach_face()
    fs.x[] = 1.;
#if TREE
  cs.refine = cs.prolongation = embed_fraction_refine;
  foreach_dimension()
    fs.x.prolongation = embed_face_fraction_refine_x;
#endif
  boundary ({cs, fs});
  restriction ({cs, fs});

  // embedded boundaries cannot be combined with (another) metric yet
  assert (is_constant (cm) || cm.i == cs.i);
  
  cm = cs;
  fm = fs;
}

References

[johansen1998]

Hans Johansen and Phillip Colella. A cartesian grid embedded boundary method for poisson’s equation on irregular domains. Journal of Computational Physics, 147(1):60-85, 1998. [ .pdf ]

See also

Usage

Tests