dr.h

Go to the documentation of this file.

/** @file dr.h
 */
/**
# Boussinesq buoyancy
 
The momentum equation of the [multilayer solver](hydro.h) becomes
\f[
\begin{aligned}
  \partial_t \left( h \mathbf{u} \right)_k + \mathbf{{\nabla}} \cdot \left(
  h \mathbf{u}  \mathbf{u} \right)_k & = - gh_k  \mathbf{{\nabla}} (\eta)
  {\color{blue} - \mathbf{{\nabla}} (h q)_k + \left[ q 
  \mathbf{{\nabla}} z \right]_k}
\end{aligned}
\f]
where the terms in blue have been added and \f$q\f$ is the (hydrostatic)
pressure deviation due to (small) density variations.
 
The density variations are described by a scalar field \f$T\f$, which is
advected by the flow, and an associated "equation of state"
\f$\Delta\rho(T)\f$ defined by the user. This gives the additional
equations
\f[
\begin{aligned}
\partial_t \left( h T \right)_k + \mathbf{{\nabla}} \cdot \left(
h \mathbf{u}  T \right)_k & = 0,\\
q(z) & = \int_0^z g \Delta \rho(T) dz
\end{aligned}
\f]
*/
 
scalar T;
 
/** @brief Event: defaults (i = 0) */
void event_defaults (void) 
{
  T = {0} /* new scalar */[nl];
  tracers = list_append (tracers, T);
}
 
/** @brief Event: cleanup (t = end, last) */
void event_cleanup (void) 
{
  delete ({T});
}
 
/**
The pressure gradient terms in blue are added to the acceleration of
the [multilayer solver](hydro.h). */
 
/** @brief Event: acceleration (i++) */
void event_acceleration (void) 
{
 
  /**
  The hydrostatic pressure deviation \f$q\f$ is stored on the interfaces
  between layers (consistently with the Keller box scheme
  discretisation, see [Popinet,
  2020](/Bibliography#popinet2020)). This gives the following vertical
  discrete integration scheme. */
 
  scalar q = {0} /* new scalar */[nl];
  
  for (int _i = 0; _i < _N; _i++) /* foreach */ {
    double ph = 0.;
    for (point.l = nl - 1; point.l >= 0; point.l--) {
      ph += G*drho(T[])*h[];
      q[] = ph;
    }
  }
 
  /**
  Once the pressure deviation is known, the terms in blue above are
  added to the face acceleration field `ha`, using the
  [pressure-gradient macro](hydro.h#horizontal-pressure-gradient). */
  
  for (int _i = 0; _i < _N; _i++) /* foreach_face */
    hpg (pg, q, 0,
   ha.x[] += pg);
 
  delete ({q});
}
 
/**
## See also
 
* [Boussinesq buoyancy for isopycnal layers](isopycnal.h)
*/