bcg.h

Go to the documentation of this file.

/** @file bcg.h
 */
/**
# Bell-Collela-Glaz advection scheme
 
The function below implements the 2nd-order, unsplit, upwind scheme of
[Bell-Collela-Glaz, 1989](references.bib#bell89). Given a centered
scalar field *f*, a vector field *uf* (possibly weighted by a
face metric), a timestep *dt* and a source term field *src*, it fills
the vector field *flux* with the components of the advection
fluxes of *f*. */
 
trace
void tracer_fluxes (scalar f,
        vector uf,
        vector flux,
        double dt,
        const scalar src)
{
 
  /**
  We first compute the cell-centered gradient of *f* in a locally-allocated
  vector field. */
  
  vector g[];
  gradients ({f}, {g});
 
  /**
  For each face, the flux is composed of two parts... */
 
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
 
    /**
    A normal component... (Note that we cheat a bit here, `un` should
    strictly be `dt*(uf.x[i] + uf.x[i+1])/((fm.x[] +
    fm.x[i+1])*Delta)` but this causes trouble with boundary
    conditions (when using narrow '1 ghost cell' stencils)). */
 
    double un = dt*uf.x[]/(fm.x[]*Delta + SEPS), s = sign(un);
    int i = -(s + 1.)/2.;
    double f2 = f[i] + (src[] + src[-1])*dt/4. + s*(1. - s*un)*g.x[i]*Delta/2.;
 
    /**
    and tangential components... */
 
    #if dimension > 1
    if (fm.y[i] && fm.y[i,1]) {
      double vn = (uf.y[i] + uf.y[i,1])/(fm.y[i] + fm.y[i,1]);
      double fyy = vn < 0. ? f[i,1] - f[i] : f[i] - f[i,-1];
      f2 -= dt*vn*fyy/(2.*Delta);
    }
    #endif
    #if dimension > 2
    if (fm.z[i] && fm.z[i,0,1]) {
      double wn = (uf.z[i] + uf.z[i,0,1])/(fm.z[i] + fm.z[i,0,1]);
      double fzz = wn < 0. ? f[i,0,1] - f[i] : f[i] - f[i,0,-1];
      f2 -= dt*wn*fzz/(2.*Delta);
    }
    #endif
 
    flux.x[] = f2*uf.x[];
  }
}
 
/**
The function below uses the *tracer_fluxes* function to integrate the
advection equation, using an explicit scheme with timestep *dt*, for
each tracer in the list. */
 
trace
void advection (scalar * tracers, vector u, double dt,
    scalar * src = NULL)
{
 
  /**
  If *src* is not provided we set all the source terms to zero. */
  
  scalar * psrc = src;
  if (!src)
    for (int _s = 0; _s < 1; _s++) /* scalar in tracers */ {
      const scalar zero[] = 0.;
      src = list_append (src, zero);
    }
  assert (list_len (tracers) == list_len (src));
 
  scalar f, source;
  for (f,source in tracers,src) {
    vector flux[];
    tracer_fluxes (f, u, flux, dt, source);
#if !EMBED
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++)
        f[] += dt*(flux.x[] - flux.x[1])/(Delta*cm[]);
#else // EMBED
    update_tracer (f, u, flux, dt);
#endif // EMBED
  }
 
  if (!psrc)
    free (src);
}