bcg.h
📄 View in API Reference (Doxygen) · View on basilisk.fr
Test cases (1): taylor-green
Bell-Collela-Glaz advection scheme
The function below implements the 2nd-order, unsplit, upwind scheme of Bell-Collela-Glaz, 1989. Given a centered scalar field *f*, a face vector field *uf* (possibly weighted by a face metric), a timestep *dt* and a source term field *src*, it fills the face vector field *flux* with the components of the advection fluxes of *f*.
trace
void tracer_fluxes (scalar f,
face vector uf,
face 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...
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, face 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 (scalar s 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) {
face vector flux[];
tracer_fluxes (f, u, flux, dt, source);
#if !EMBED
foreach()
foreach_dimension()
f[] += dt*(flux.x[] - flux.x[1])/(Delta*cm[]);
#else // EMBED
update_tracer (f, u, flux, dt);
#endif // EMBED
}
if (!psrc)
free (src);
}