mytracer.h

Tracer advection event

Tracer advection event

This event integrates advection equations of the form \partial_tf_i+\mathbf{u_f}\cdot\nabla f_i=0 where \mathbf{u_f} is the velocity field and f_i are a list of passive tracers.

The tracers list is defined elsewhere (typically by the user), the face vector field uf and the timestep dt are defined by a solver.

extern scalar * tracers;
extern face vector uf;
extern double dt;

extern scalar cs, csm1;
extern bool emerged;

On adaptive meshes, tracers need to use linear interpolation (rather than the default bilinear interpolation) to ensure conservation when refining cells.

#if TREE
event defaults (i = 0) {
  for (scalar s in tracers) {
#if EMBED
    s.refine = s.prolongation = refine_embed_linear;
#else
    s.refine  = refine_linear;
#endif
    s.restriction = restriction_volume_average;
  }
}
#endif

The integration is performed using the Bell-Collela-Glaz scheme.

#include "bcg.h"

event tracer_advection (i++,last) {
  advection (tracers, uf, dt);
}

Diffusion can be added by overloading this hook.

event tracer_diffusion (i++,last);

We update the tracer in the solid and emerged cells, using the boundary condition at the embedded boundary.

event advection_term (i++,last)
{

We first make sure not to use any values in newly emerged cells by setting the flag emerged to false.

  emerged = false;
  boundary (tracers);

In the solid cells, we set the scalar field s to 0.

  for (scalar s in tracers)
    foreach()
      if (cs[] <= 0.) 
	s[] = 0.;

In the emerged cells, we use the solid boundary condition to update the scalar field.

  for (scalar s in tracers)
    foreach() {
      if (csm1[] <= 0. && cs[] > 0.) {

	// Normal emerged cell
	if (cs[] < 1.) {

	  // Cell centroid, barycenter and normal of the embedded fragment
	  coord c, b, n;
	  embed_geometry (point, &b, &n);
	  double alpha = plane_alpha (cs[], n);
	  plane_center (n, alpha, cs[], &c); // Different from line_area_center

	  // Boundary condition on the embedded boundary
	  bool dirichlet = true;
	  double sb = (s.boundary[embed] (point, point, s, &dirichlet));
	  if (!dirichlet) {
	    double coef = 0.;
	    sb = neumann_scalar (point, s, cs, n, b, sb, &coef);
	    // s[] is undefined here so we use the average over the neighbors
	    int navg = 0;
	    double savg = 0.;
	    foreach_neighbor(1)
	      if (cs[] > 0. && (emerged || csm1[] > 0.)) {
		navg += 1;
		savg += s[];
	      }
	    sb += coef*(savg/(navg + SEPS));
	  }

	  // Emerged cell value
	  s[] = embed_extrapolate (point, s, cs, n, c, sb);
	}

	// Pathological emerged cell (cs = 1)
	else {
	  int navg = 0;
	  double savg = 0.;
	  foreach_neighbor(1)
	    if (cs[] > 0. && (emerged || csm1[] > 0.)) {
	      navg += 1;
	      savg += s[];
	    }
	  s[] = savg/(navg + SEPS);
	}
      }
    }

Before using the boundary function, we set the emerged flag to true to indicate that all emerged cells have been updated and can now be used.

  emerged = true;
  boundary (tracers);
}