centered.h

Go to the documentation of this file.

/** @file centered.h
 */
/**
# Incompressible Navier--Stokes solver (centered formulation)
 
We wish to approximate numerically the incompressible,
variable-density Navier--Stokes equations
\f[
\partial_t\mathbf{u}+\nabla\cdot(\mathbf{u}\otimes\mathbf{u}) = 
\frac{1}{\rho}\left[-\nabla p + \nabla\cdot(2\mu\mathbf{D})\right] + 
\mathbf{a}
\f]
\f[
\nabla\cdot\mathbf{u} = 0
\f]
with the deformation tensor 
\f$\mathbf{D}=[\nabla\mathbf{u} + (\nabla\mathbf{u})^T]/2\f$.
 
The scheme implemented here is close to that used in Gerris ([Popinet,
2003](/src/references.bib#popinet2003), [Popinet,
2009](/src/references.bib#popinet2009), [Lagrée et al,
2011](/src/references.bib#lagree2011)).
 
We will use the generic time loop, a CFL-limited timestep, the
Bell-Collela-Glaz advection scheme and the implicit viscosity
solver. If embedded boundaries are used, a different scheme is used
for viscosity. */
 
#include "run.h"
#include "timestep.h"
#include "bcg.h"
#if EMBED
# include "viscosity-embed.h"
#else
# include "viscosity.h"
#endif
 
/**
The primary variables are the centered pressure field \f$p\f$ and the
centered velocity field \f$\mathbf{u}\f$. The centered vector field
\f$\mathbf{g}\f$ will contain pressure gradients and acceleration terms.
 
We will also need an auxilliary face velocity field \f$\mathbf{u}_f\f$ and
the associated centered pressure field \f$p_f\f$. */
 
scalar p[];
vector u[], g[];
scalar pf[];
vector uf[];
 
/**
In the case of variable density, the user will need to define both the
face and centered specific volume fields (\f$\alpha\f$ and \f$\alpha_c\f$
respectively) i.e. \f$1/\rho\f$. If not specified by the user, these
fields are set to one i.e. the density is unity.
 
Viscosity is set by defining the face dynamic viscosity \f$\mu\f$; default
is zero.
 
The face field \f$\mathbf{a}\f$ defines the acceleration term; default is
zero.
 
The statistics for the (multigrid) solution of the pressure Poisson
problems and implicit viscosity are stored in *mgp*, *mgpf*, *mgu*
respectively. 
 
If *stokes* is set to *true*, the velocity advection term
\f$\nabla\cdot(\mathbf{u}\otimes\mathbf{u})\f$ is omitted. This is a
reference to [Stokes flows](http://en.wikipedia.org/wiki/Stokes_flow)
for which inertia is negligible compared to viscosity. */
 
const vector mu = zerof, a = zerof, alpha = unityf;
const scalar rho = unity;
mgstats mgp = {0}, mgpf = {0}, mgu = {0};
bool stokes = false;
 
/**
## Boundary conditions
 
For the default symmetric boundary conditions, we need to ensure that
the normal component of the velocity is zero after projection. This
means that, at the boundary, the acceleration \f$\mathbf{a}\f$ must be
balanced by the pressure gradient. Taking care of boundary orientation
and staggering of \f$\mathbf{a}\f$, this can be written */
 
#if EMBED
# define neumann_pressure(i) (alpha.n[i] ? a.n[i]*fm.n[i]/alpha.n[i] :  \
            a.n[i]*rho[]/(cm[] + SEPS))
#else
# define neumann_pressure(i) (a.n[i]*fm.n[i]/alpha.n[i])
#endif
 
p[right] = neumann (neumann_pressure(ghost));
p[left]  = neumann (- neumann_pressure(0));
 
#if AXI
uf.n[bottom] = 0.;
uf.t[bottom] = dirichlet(0); // since uf is multiplied by the metric which
                             // is zero on the axis of symmetry
p[top]    = neumann (neumann_pressure(ghost));
#else // !AXI
#  if dimension > 1
p[top]    = neumann (neumann_pressure(ghost));
p[bottom] = neumann (- neumann_pressure(0));
#  endif
#  if dimension > 2
p[front]  = neumann (neumann_pressure(ghost));
p[back]   = neumann (- neumann_pressure(0));
#  endif
#endif // !AXI
 
/**
For [embedded boundaries on trees](/src/embed-tree.h), we need to
define the pressure gradient for prolongation of pressure close to
embedded boundaries. */
 
#if TREE && EMBED
void pressure_embed_gradient (Point point, scalar p, coord * g)
{
  for (int _d = 0; _d < dimension; _d++)
    g->x = rho[]/(cm[] + SEPS)*(a.x[] + a.x[1])/2.;
}
#endif // TREE && EMBED
 
/**
## Initial conditions */
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) 
{
 
  /**
  We reset the multigrid parameters to their default values. */
  
  mgp = (mgstats){0};
  mgpf = (mgstats){0};
  mgu = (mgstats){0};  
  
  CFL = 0.8;
 
  /**
  The pressures are never dumped. */
 
  p.nodump = pf.nodump = true;
  
  /**
  The default density field is set to unity (times the metric and the
  solid factors). */
 
  if (alpha.x.i == unityf.x.i) {
    alpha = fm;
    rho = cm;
  }
  else if (!is_constant(alpha.x)) {
    vector alphav = alpha;
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      alphav.x[] = fm.x[];
  }
 
  /**
  On trees, refinement of the face-centered velocity field needs to
  preserve the divergence-free condition. */
 
#if TREE
  uf.x.refine = refine_face_solenoidal;
 
  /**
  When using [embedded boundaries](/src/embed.h), the restriction and
  prolongation operators need to take the boundary into account. */
 
#if EMBED
  uf.x.refine = refine_face;
  for (int _d = 0; _d < dimension; _d++)
    uf.x.prolongation = refine_embed_face_x;
  for (int _s = 0; _s < 1; _s++) /* scalar in {p, pf, u, g} */ {
    s.restriction = restriction_embed_linear;
    s.refine = s.prolongation = refine_embed_linear;
    s.depends = list_add (s.depends, cs);
  }
  for (int _s = 0; _s < 1; _s++) /* scalar in {p, pf} */
    s.embed_gradient = pressure_embed_gradient;
#endif // EMBED
#endif // TREE
 
  /**
  We set the dimensions of the velocity field. */
 
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++)
      /* dim: u.x[] == Delta/t */;
}
 
 
/**
We had some objects to display by default. */
 
/** @brief Event: default_display (i = 0) */
void event_default_display (void) 
  display ("squares (color = 'u.x', spread = -1);");
 
/**
After user initialisation, we initialise the face velocity and fluid
properties. */
 
double dtmax;
 
/** @brief Event: init (i = 0) */
void event_init (void) 
{
  /* trash: uf */;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    uf.x[] = fm.x[]*face_value (u.x, 0);
 
  /**
  We update fluid properties. */
 
  event ("properties");
 
  /**
  We set the initial timestep (this is useful only when restoring from
  a previous run). */
 
  dtmax = DT;
  event ("stability");
}
 
/**
## Time integration
 
The timestep for this iteration is controlled by the CFL condition,
applied to the face centered velocity field \f$\mathbf{u}_f\f$; and the
timing of upcoming events. */
 
event set_dtmax (i++,last) dtmax = DT;
 
/** @brief Event: stability (i++,last) */
void event_stability (void) {
  dt = dtnext (stokes ? dtmax : timestep (uf, dtmax));
}
 
/**
If we are using VOF or diffuse tracers, we need to advance them (to
time \f$t+\Delta t/2\f$) here. Note that this assumes that tracer fields
are defined at time \f$t-\Delta t/2\f$ i.e. are lagging the
velocity/pressure fields by half a timestep. */
 
event vof (i++,last);
event tracer_advection (i++,last);
event tracer_diffusion (i++,last);
 
/**
The fluid properties such as specific volume (fields \f$\alpha\f$ and
\f$\alpha_c\f$) or dynamic viscosity (face field \f$\mu_f\f$) -- at time
\f$t+\Delta t/2\f$ -- can be defined by overloading this event. */
 
event properties (i++,last);
 
/**
### Predicted face velocity field
 
For second-order in time integration of the velocity advection term
\f$\nabla\cdot(\mathbf{u}\otimes\mathbf{u})\f$, we need to define the face
velocity field \f$\mathbf{u}_f\f$ at time \f$t+\Delta t/2\f$. We use a version
of the Bell-Collela-Glaz [advection scheme](/src/bcg.h) and the
pressure gradient and acceleration terms at time \f$t\f$ (stored in vector
\f$\mathbf{g}\f$). */
 
void prediction()
{
  vector du;
  for (int _d = 0; _d < dimension; _d++) {
    scalar s = {0} /* new scalar */;
    du.x = s;
  }
 
  if (u.x.gradient)
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++) {
#if EMBED
        if (!fs.x[] || !fs.x[1])
    du.x[] = 0.;
  else
#endif
    du.x[] = u.x.gradient (u.x[-1], u.x[], u.x[1])/Delta;
      }
  else
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++) {
#if EMBED
        if (!fs.x[] || !fs.x[1])
    du.x[] = 0.;
  else
#endif
    du.x[] = (u.x[1] - u.x[-1])/(2.*Delta);
    }
 
  /* trash: uf */;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
    double un = dt*(u.x[] + u.x[-1])/(2.*Delta), s = sign(un);
    int i = -(s + 1.)/2.;
    uf.x[] = u.x[i] + (g.x[] + g.x[-1])*dt/4. + s*(1. - s*un)*du.x[i]*Delta/2.;
    #if dimension > 1
    if (fm.y[i,0] && fm.y[i,1]) {
      double fyy = u.y[i] < 0. ? u.x[i,1] - u.x[i] : u.x[i] - u.x[i,-1];
      uf.x[] -= dt*u.y[i]*fyy/(2.*Delta);
    }
    #endif
    #if dimension > 2
    if (fm.z[i,0,0] && fm.z[i,0,1]) {
      double fzz = u.z[i] < 0. ? u.x[i,0,1] - u.x[i] : u.x[i] - u.x[i,0,-1];
      uf.x[] -= dt*u.z[i]*fzz/(2.*Delta);
    }
    #endif
    uf.x[] *= fm.x[];
  }
 
  delete ((scalar *){du});
}
 
/**
### Advection term
 
We predict the face velocity field \f$\mathbf{u}_f\f$ at time \f$t+\Delta
t/2\f$ then project it to make it divergence-free. We can then use it to
compute the velocity advection term, using the standard
Bell-Collela-Glaz advection scheme for each component of the velocity
field. */
 
/** @brief Event: advection_term (i++,last) */
void event_advection_term (void) 
{
  if (!stokes) {
    prediction();
    mgpf = project (uf, pf, alpha, dt/2., mgpf.nrelax);
    advection ((scalar *){u}, uf, dt, (scalar *){g});
  }
}
 
/**
### Viscous term
 
We first define a function which adds the pressure gradient and
acceleration terms. */
 
static void correction (double dt)
{
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++)
      u.x[] += dt*g.x[];
}
 
/**
The viscous term is computed implicitly. We first add the pressure
gradient and acceleration terms, as computed at time \f$t\f$, then call
the implicit viscosity solver. We then remove the acceleration and
pressure gradient terms as they will be replaced by their values at
time \f$t+\Delta t\f$. */
 
/** @brief Event: viscous_term (i++,last) */
void event_viscous_term (void) 
{
  if (constant(mu.x) != 0.) {
    correction (dt);
    mgu = viscosity (u, mu, rho, dt, mgu.nrelax);
    correction (-dt);
  }
 
  /**
  We reset the acceleration field (if it is not a constant). */
 
  if (!is_constant(a.x)) {
    vector af = a;
    /* trash: af */;
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      af.x[] = 0.;
  }
}
 
/**
### Acceleration term
 
The acceleration term \f$\mathbf{a}\f$ needs careful treatment as many
equilibrium solutions depend on exact balance between the acceleration
term and the pressure gradient: for example Laplace's balance for
surface tension or hydrostatic pressure in the presence of gravity.
 
To ensure a consistent discretisation, the acceleration term is
defined on faces as are pressure gradients and the centered combined
acceleration and pressure gradient term \f$\mathbf{g}\f$ is obtained by
averaging. 
 
The (provisionary) face velocity field at time \f$t+\Delta t\f$ is
obtained by interpolation from the centered velocity field. The
acceleration term is added. */
 
/** @brief Event: acceleration (i++,last) */
void event_acceleration (void) 
{
  /* trash: uf */;
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    uf.x[] = fm.x[]*(face_value (u.x, 0) + dt*a.x[]);
}
 
/**
## Approximate projection
 
This function constructs the centered pressure gradient and
acceleration field *g* using the face-centered acceleration field *a*
and the cell-centered pressure field *p*. */
 
void centered_gradient (scalar p, vector g)
{
 
  /**
  We first compute a face field \f$\mathbf{g}_f\f$ combining both
  acceleration and pressure gradient. */
 
  vector gf[];
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    gf.x[] = fm.x[]*a.x[] - alpha.x[]*(p[] - p[-1])/Delta;
 
  /**
  We average these face values to obtain the centered, combined
  acceleration and pressure gradient field. */
 
  /* trash: g */;
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++)
      g.x[] = (gf.x[] + gf.x[1])/(fm.x[] + fm.x[1] + SEPS);
}
 
/**
To get the pressure field at time \f$t + \Delta t\f$ we project the face
velocity field (which will also be used for tracer advection at the
next timestep). Then compute the centered gradient field *g*. */
 
/** @brief Event: projection (i++,last) */
void event_projection (void) 
{
  mgp = project (uf, p, alpha, dt, mgp.nrelax);
  centered_gradient (p, g);
 
  /**
  We add the gradient field *g* to the centered velocity field. */
 
  correction (dt);
}
 
/**
Some derived solvers need to hook themselves at the end of the
timestep. */
 
event end_timestep (i++, last);
 
/**
## Adaptivity
 
After mesh adaptation fluid properties need to be updated. When using
[embedded boundaries](/src/embed.h) the fluid fractions and face
fluxes need to be checked for inconsistencies. */
 
#if TREE
/** @brief Event: adapt (i++,last) */
void event_adapt (void) {
#if EMBED
  fractions_cleanup (cs, fs);
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    if (uf.x[] && !fs.x[])
      uf.x[] = 0.;
#endif
  event ("properties");
}
#endif
 
/**
## See also
 
* [Double projection](double-projection.h)
* [Performance monitoring](perfs.h)
*/