all-mach.h

Go to the documentation of this file.

/** @file all-mach.h
 */
/**
# An "all Mach" flow solver
 
We wish to solve the generic momentum equation
\f[
\partial_t\mathbf{q} + \nabla\cdot(\mathbf{q}\mathbf{u}) = 
- \nabla p + \nabla\cdot(\mu\nabla\mathbf{u}) + \rho\mathbf{a}
\f]
with \f$\mathbf{q}=\rho\mathbf{u}\f$ the momentum, \f$\mathbf{u}\f$ the
velocity, \f$\rho\f$ the density, \f$\mu\f$ the dynamic viscosity, \f$p\f$ the
pressure and \f$\mathbf{a}\f$ an
acceleration. The pressure is defined through an equation of state and
verifies the evolution equation
\f[
\partial_t p + \mathbf{u}\cdot\nabla p = -\rho c^2\nabla\cdot\mathbf{u}
\f]
with \f$c\f$ the speed of sound. By default the solver sets \f$c=\infty\f$,
\f$\rho=1\f$ and the pressure equation reduces to 
\f[
\nabla\cdot\mathbf{u} = 0
\f]
 
The advection of momentum is not performed by this solver (so that
different schemes can be used) i.e. in the end, by default, we solve
the incompressible (linearised) Euler equations with a projection
method.
 
We build the solver using the generic time loop and the implicit
viscous solver (which includes the multigrid Poisson--Helmholtz
solver). */
 
#include "run.h"
#include "timestep.h"
#include "viscosity.h"
 
/**
The primitive variables are the momentum \f$\mathbf{q}\f$, pressure \f$p\f$,
density \f$\rho\f$, (face) specific volume \f$\alpha=1/\rho\f$, (face) dynamic
viscosity \f$\mu\f$ (which is zero by default) and (face) velocity field
\f$\mathbf{u}_f\f$. */
 
vector q[];
scalar p[];
vector uf[];
const vector alpha = unityf, mu = zerof;
const scalar rho = unity;
 
/**
The equation of state is defined by the pressure field *ps* and \f$\rho
c^2\f$. In the incompressible limit \f$\rho c^2\rightarrow\infty\f$. Rather
than trying to compute this limit, we set both fields to zero and
check for this particular case when computing the pressure (see
below). This means that by default the fluid is incompressible. */
 
scalar ps[];
const scalar rhoc2 = zeroc;
const vector a = zerof;
 
/**
We store the combined pressure gradient and acceleration field in
*g*. */
 
vector g[];
mgstats mgp = {0}, mgu = {0};
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) {
 
  /**
  The default density field is set to unity (times the metric). */
 
  if (alpha.x.i == unityf.x.i)
    alpha = fm;
 
  if (rho.i == unity.i)
    rho = cm;
  
  /**
  We reset the multigrid parameters to their default values. */
  
  mgp = (mgstats){0};
  mgu = (mgstats){0};    
}
 
/** @brief Event: init (i = 0) */
void event_init (void) {
 
  /**
  The face velocity field is obtained by simple linear interpolation
  from the momentum field. We make sure that the specific volume
  \f$\alpha\f$ is defined by calling the "properties" event (see
  below). */
  
  event ("properties");
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    uf.x[] = alpha.x[]*(q.x[] + q.x[-1])/2.;
}
 
/**
The timestep is computed using the CFL condition on the face velocity
field. */
 
double dtmax;
 
event set_dtmax (i++,last) dtmax = DT;
 
/** @brief Event: stability (i++,last) */
void event_stability (void) {
  dt = dtnext (timestep (uf, dtmax));
}
 
/**
Tracers (including momentum \f$\mathbf{q}\f$) are advected by these events. */
 
event vof (i++,last);
event tracer_advection (i++,last);
 
/**
The equation of state (i.e. fields \f$\alpha\f$, \f$\rho\f$, \f$\rho c^2\f$ and
*ps*) is defined by this event. */
 
/** @brief Event: properties (i++,last) */
void event_properties (void) 
{
  
  /**
  If the acceleration is not constant, we reset it to zero. */
  
  if (!is_constant(a.x)) {
    vector af = a;
    for (int _i = 0; _i < _N; _i++) /* foreach_face */
      af.x[] = 0.;
  }
}
 
/**
This event can be overloaded to add acceleration terms. */
 
event acceleration (i++, last);
 
/**
The equation for the pressure is a Poisson--Helmoltz problem which we
will solve with the [multigrid solver](poisson.h). The statistics for
the solver will be stored in *mgp* (resp. *mgu* for the viscosity
solver). */
 
/** @brief Event: pressure (i++, last) */
void event_pressure (void) 
{
 
  /**
  If the viscosity is not zero, we use the implicit viscosity solver
  to obtain the velocity field at time \f$t + \Delta t\f$. The pressure
  term is taken into account using the pressure gradient at time
  \f$t\f$. A provisionary momentum (without the pressure gradient) is then
  built from the velocity field. */
  
  if (constant(mu.x) != 0.) {
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++)
        q.x[] = (q.x[] + dt*g.x[])*cm[]/rho[];
    mgu = viscosity (q, mu, rho, dt, mgu.nrelax);
    for (int _i = 0; _i < _N; _i++) /* foreach */
      for (int _d = 0; _d < dimension; _d++)
        q.x[] = q.x[]*rho[]/cm[] - dt*g.x[];
  }
 
  /**
  We first define a temporary face velocity field \f$\mathbf{u}_\star\f$
  using simple averaging from \f$\mathbf{q}_{\star}\f$, \f$\alpha_{n+1}\f$ and
  the acceleration term. */
 
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    uf.x[] = alpha.x[]*(q.x[] + q.x[-1])/2. + dt*fm.x[]*a.x[];
 
  /**
  The evolution equation for the pressure is
  \f[\partial_tp +\mathbf{u} \cdot \nabla p = - \rho c^2 \nabla \cdot \mathbf{u}\f]
  with \f$\rho\f$ the density and \f$c\f$ the speed of sound. Following the
  classical [projection
  method](navier-stokes/centered.h#approximate-projection) for
  incompressible flows, we set
  \f[
  \mathbf{u}_{n + 1} = \mathbf{u}_{\star} - \Delta t (\alpha\nabla p)_{n+1}
  \f]
  The evolution equation for the pressure can then be discretised as
  \f[
  \frac{p_{n + 1} - p_n}{\Delta t} +\mathbf{u}_n \cdot \nabla p_n = 
     - \rho c^2_{n + 1} \nabla \cdot \mathbf{u}_{n + 1}
  \f]
  which gives, after some manipulations, the Poisson--Helmholtz equation
  \f[
  \lambda_{n + 1} p_{n + 1} + \nabla \cdot \left( \alpha \nabla p
  \right)_{n + 1} = \lambda_{n + 1} p_{\star} + \frac{1}{\Delta t} \nabla \cdot
  \mathbf{u}_{\star}
  \f]
  with
  \f[
  p_{\star} = p_n - \Delta t\mathbf{u}_n \cdot \nabla p_n
  \f]
  and
  \f[
  \lambda = \frac{- 1}{\Delta t^2 \rho c^2}
  \f]
  */
  
  scalar lambda = rhoc2, rhs = ps;
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
 
    /**
    We compute \f$\lambda\f$ and the corresponding term in the
    right-hand-side of the Poisson--Helmholtz equation. */
 
    if (constant(lambda) == 0.)
      rhs[] = 0.;
    else {
      lambda[] = - cm[]/(sq(dt)*rhoc2[]);
      rhs[] = lambda[]*ps[];
    }
      
    /**
    We add the divergence of the velocity field to the right-hand-side. */
 
    double div = 0.;
    for (int _d = 0; _d < dimension; _d++)
      div += uf.x[1] - uf.x[];
    rhs[] += div/(dt*Delta);
  }
  
  /**
  The Poisson--Helmholtz solver is called with a [definition of the
  tolerance](poisson.h#project) identical to that used for
  incompressible flows. */
  
  mgp = poisson (p, rhs, alpha, lambda, tolerance = TOLERANCE/sq(dt));
  
  /**
  The pressure gradient is applied to \f$\mathbf{u}_\star\f$ to obtain the
  face velocity field at time \f$n + 1\f$. 
  
  We also compute the combined face pressure gradient and acceleration
  field *gf*. */
 
  vector gf[];
  for (int _i = 0; _i < _N; _i++) /* foreach_face */ {
    double dp = alpha.x[]*(p[] - p[-1])/Delta;
    uf.x[] -= dt*dp;
    gf.x[] = a.x[] - dp/fm.x[];
  }
 
  /**
  And finally we apply the pressure gradient/acceleration term to the
  flux/momentum. We also store the centered, combined pressure
  gradient and acceleration field *g*. */
  
  for (int _i = 0; _i < _N; _i++) /* foreach */
    for (int _d = 0; _d < dimension; _d++) {
      g.x[] = rho[]*(gf.x[] + gf.x[1])/(2.*cm[]);
      q.x[] += dt*g.x[];
    }
}
 
/**
Some derived solvers need to hook themselves at the end of the
timestep. */
 
event end_timestep (i++, last);
 
/**
After mesh adaptation fluid properties need to be updated. */
 
#if TREE
/** @brief Event: adapt (i++,last) */
void event_adapt (void) {
  event ("properties");
}
#endif