%-------------------------------------------------------------------------------

% This file is part of Code_Saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2021 EDF S.A.
%
% This program is free software; you can redistribute it and/or modify it under
% the terms of the GNU General Public License as published by the Free Software
% Foundation; either version 2 of the License, or (at your option) any later
% version.
%
% This program is distributed in the hope that it will be useful, but WITHOUT
% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
% details.
%
% You should have received a copy of the GNU General Public License along with
% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
% Street, Fifth Floor, Boston, MA 02110-1301, USA.

%-------------------------------------------------------------------------------

%-------------------------------------------------------------------------------
\section{Time discretisation of a transport equation}

At first, the physical properties of the flow are computed (density,
viscosity, specific heat \emph{etc.}): indeed, they may depend upon the variables
(such as the temperature for example).

The time scheme is a $\theta$-scheme:
\begin{equation}
\left\{%
\begin{array}{ll}
\theta = 1 & \text{for an implicit first order Euler scheme}, \\
\theta = \frac{1}{2} & \text{for second order Crank-Nicolson scheme}.
\end{array}
\right.
\end{equation}

For the second order scheme, the time step is assumed to be constant.

If required, the equations for the turbulent variables are solved (turbulent
kinetic energy and dissipation or Reynolds stresses and dissipation), using
a $\theta$-scheme again. For the $k-\varepsilon$ model, an additional step
is carried out to couple the source terms. For the Reynolds stress model,
the variables (turbulent stresses and dissipation) are solved sequentially,
without coupling.

Next, the equations for the \emph{scalars} (enthalpy, temperature, tracers,
concentrations, mass fractions...) are solved, also with a $\theta$-scheme.

Finally, all the variables are updated and another time step may start.

The general equation for advection (valid for the velocity components, the
turbulent variables and the scalars) is re-written as follows in a condensed
form; the mass equation ($\frac{\partial \rho } {\partial t}+ \dive(\rho
\vect{u}) = \Gamma$ see Equation \eqref{eq:goveqn:mass}) has been used to split the time derivative:
%
\begin{equation}\label{eq:timstp:transport_eq_scalar}
\rho \der{\varia}{t} + \grad \varia \cdot \left(\rho \vect{u} \right)
- \dive \left(K \grad \varia \right) = S_{i}\left(\Phi, \, \varphi\right) \varia + S_{e}\left(\Phi,\varphi \right)
+ \Gamma \left(\varia^{in} - \varia \right).
\end{equation}
In Equation \eqref{eq:timstp:transport_eq_scalar}, $\Phi$ represents the physical properties such as $(\rho,K,\mu_{t},...)$,
$\varphi$ represents the variables of the problem such as $(\vect{u},k,\epsilon,...)$,
$S_{i}\left(\Phi, \, \varphi \right) \varia$ is the linear part of the source terms
and
$S_{e} \left(\Phi, \, \varphi \right)$  includes all other source terms.

Therefore, four different time steppings are used, they all define the time at which the quantities are evaluated:
%
\begin{enumerate}[ label=\roman{*}/, ref=(\roman{*})]
\item $\theta$ is the time stepping applied to the variable $\varia^{n+\theta} \equiv \theta \varia^{n+1} + \left( 1 - \theta \right) \varia^n$,
\item $\theta_\Phi$ is the time stepping applied to the physical properties,
\item $\theta_F$ is the time stepping applied to the mass flux,
\item $\theta_S$ is the time stepping applied to the source terms,
\end{enumerate}
%

If $\theta=1/2$, or if an extrapolation is used, the time step $\Delta t$ is
constant in time and uniform in space.

\subsection{Physical properties $\Phi$}

The physical properties of the flow (density, viscosity, specific heat...)
are:

\begin{itemize}
\item either explicit, defined at the time step $n$.
\item or extrapolated at $n+\theta _{\Phi }$ using the Adam-Bashforth
time scheme (in this case, the time step is assumed to be constant).
\end{itemize}

Under a more general form, this reads:
\begin{equation}
 \Phi^{n+\theta_{\Phi}} \equiv(1+\theta_{\Phi})\,\Phi^{n}- \theta_{\Phi}\,\Phi^{n-1},
\end{equation}

\begin{equation}
\left\{%
\begin{array}{ll}
\theta_{\Phi} = 0 & \text{standard explicit formulation} ,\\
\theta_{\Phi} = 1/2 & \text{second order extrapolation at } n+1/2, \\
\theta_{\Phi} = 1 & \text{first order extrapolation at } n+1.
\end{array}
\right.
\end{equation}

\subsection{Mass flux}

\hypertarget{massflux}{}

For the mass flux, three time schemes are available. The mass flux may be:
%
\begin{itemize}
\item explicit, taken at time step $n$ for the momentum equations and
updated with its value at time step $n+1$ for the equations for turbulence
and scalars (standard scheme).

\item explicit, taken at time step $n$ for the momentum equations and
also for the equations for turbulence and scalars.

\item taken at $n+\theta_{F}$ (second order if $\theta_{F}=1/2$). To
solve the momentum equations, $\left(\rho\,\vect{u}\right)^{n-2+\theta_{F}}$ and
$\left(\rho\,\vect{u}\right)^{n-1+\theta_{F}}$ are known. Hence, the value at
$n+\theta_{F}$ is obtained as a result of the following extrapolation:
\begin{equation}
\left(\rho\vect{u}\right)^{n+\theta_{F}}
= 2 \left(\rho \vect{u}\right)^{n-1+\theta_{F}}
- \left(\rho \vect{u}\right)^{n-2+\theta_{F}}.
\end{equation}
At the end of this phase (after the pressure correction step),
$\left(\rho\,\vect{u}\right)^{n+1}$ is known and the following interpolation is used to
determine the mass flux at $n+\theta_{F}$ that will be adopted for the
equations for turbulence and scalars:
\begin{equation}
\left(\rho\,\vect{u}\right)^{n+\theta_{F}}= \dfrac{1}{2-\theta_{F}}
   \left(\rho \vect{u}\right)^{n+1}
+\frac{1-\theta_{F}}{2-\theta_{F}} \left(\rho \vect{u}\right)^{n-1+\theta_{F}}.
\end{equation}
\end{itemize}

See the \doxygenfile{schtmp_8f90.html}{programmers reference of the dedicated subroutine} for further details
on the computation of the mass flux.

\subsection{Source terms}

As for the physical properties, the \textbf{explicit} source terms are:

\begin{itemize}
\item explicit:
\begin{equation}
\left[ S_{e} \left(\Phi , \, \varphi \right) \right]^{n}=S_{e}\left(\Phi ^{n+\theta_{\Phi }}, \, \varphi^{n} \right),
\end{equation}

\item extrapolated at $n+\theta _{S}$ using the Adams-Bashforth scheme:
\begin{equation}
\left[ S_{e} \left(\Phi , \, \varphi \right) \right]^{n+\theta _{S}}=\left(1+\theta _{S}\right)\,
S_{e}\left(\Phi^{n}, \, \varphi^{n}\right)-\theta _{S}\,S_{e} \left(\Phi ^{n-1} , \,\varphi ^{n-1}\right) .
\end{equation}
\end{itemize}

By default, to be consistent and preserve the order of convergence in time,
the implicit source terms are discretized with the same scheme as that is used
for convection-diffusion of the unknown considered, \emph{i.e.} taken at
$n+\theta $:
\begin{equation}
\left[ S_{i} \left(\Phi , \, \varphi \right) \varia \right]^{n+\theta }
= S_{i}(\Phi^{n+\theta _{\Phi}} , \, \varphi^{n}) \,\left[ \theta \varia^{n+1}+\left(1-\theta \right)\varia^{n} \right].
\end{equation}

\begin{remark}
The \textbf{implicit} source terms taken also at $n+\theta $ for $\theta
_{S}\neq 0$, while for $\theta _{S}=0$, the implicit source terms are taken
at $n+1$ , this to enhance stability.
\end{remark}

\subsection{General time discretized form}

For the sake of clarity, it is assumed hereafter that, unless otherwise
explicitly stated, the mass flux is taken at $n+\theta_F$ and the physical
properties are taken at $n+\theta_\Phi$, with $\theta_F$ and $\theta_\Phi$
dependent upon the specific schemes selected for the mass flux and the
physical properties respectively and all $\theta$s are denoted by $\theta$.

Under a general form, the discrete counterpart of Equation \eqref{eq:timstp:transport_eq_scalar} at
$n+\theta$ reads:
\begin{equation}
\displaystyle \dfrac{\rho}{\Delta t} \left(\varia^{n+1}-\varia^{n}\right)
+ \grad\varia^{n+\theta} \cdot \left(\rho \vect{u} \right)
- \dive \left( K \grad \varia^{n+\theta} \right) =
\left[ S_{i}(\Phi, \, \varphi) \varia \right]^{n+\theta} + \left[S_{e}(\Phi, \, \varphi) \right]^{n+\theta} .
\end{equation}

Using the standard $\theta$-scheme $\varia^{n+\theta}=\theta \varia^{n+1}+\left(1-\theta \right)
\varia^{n}$, the equation reads:
\begin{equation}
\begin{array}{r @{\,} c @{\,} l}
\displaystyle \dfrac{\rho}{\Delta t} \left( \varia^{n+1}-\varia^{n} \right)
+ \theta \left[
\grad \varia^{n+1} \cdot \left( \rho \vect{u} \right)
- \dive \left( K \grad \varia^{n+1} \right)   \right]
&=&
- \left(1-\theta\right) \left[
\grad \varia^{n} \cdot \left( \rho \vect{u} \right)
- \dive \left( K \grad \varia^{n}\right)
\right]
 \\
&+& S_{i}\left(\Phi , \,\varphi^n\right)\left[\theta \varia^{n+1}+ \left(1-\theta \right) \varia^{n}\right]
+ \left[ S_{e}(\Phi, \, \varphi)\right]^{n+\theta}.
\end{array}
\end{equation}

For numerical reasons, the system is solved in an iterative and incremental
manner, with the help of the series $\delta \varia^{n+1}_{k+1}=\varia^{n+1}_{k+1}-\varia^{n+1}_{k}$
(with, by definition, $\varia^{n+1}_0=\varia^{n}$).
More theoretical details of such an iterative process are given in \S \ref{sec:algebrea_iterative_process}.

%-------------------------------------------------------------------------------
\section{Pressure-based velocity-pressure solver}
The aim of this section is to describe how Navier Stokes equations are solved for an incompressible
or weakly compressible (dilatable or Low Mach algorithm) combined with an
implicit Euler time stepping or a second order Crank Nicolson time stepping.

The set of equations to be solved is
\begin{equation}\label{eq:timstp:navier_stokes}
\left\{\begin{array}{l}
\displaystyle \rho \der{\vect{u}}{t} + \gradt \, \vect{u} \cdot \left( \rho\vect{u}  \right)
=\divv(\tens{\sigma}) - \divv(\rho \tens{R}) + \vect{ST}_u - \tens{K}\,\vect{u} + \rho \vect{g} + \Gamma \left( \vect{u}^{in} - \vect{u} \right), \\
\dfrac{\partial \rho}{\partial t} + \dive\left( \rho \vect{u} \right) = \Gamma ,
\end{array}\right.
\end{equation}
%
where $\rho$ is the density field, $\vect{u}$ is the velocity field to be solved,
$[\,\vect{ST}_u-\tens{K}\,\vect{u} + \rho \vect{g} \,]$  are sources terms (note that $\tens{K}$ is expected to be a positive definite tensor),
 $\tens{\sigma}$ is the stress tensor, composed of the viscous tensor $\tens{\tau}$ and of the pressure field as follows
 %
\begin{equation}
\left\{\begin{array}{r c l}
\tens{\sigma} &=& \tens{\tau} - P\tens{Id} , \\%FIXME +Rij
\tens{\tau} & =& 2 \mu  \tens{S} +   \left( \kappa - \dfrac{2}{3}\mu \right)   \trace \left( \tens{S} \right) \tens{1} ,\\
\tens{S} &=& \dfrac{1}{2} \left( \gradt \, \vect{u}+\transpose{\gradt \, \vect{u}} \right) .
\end{array}\right.
\end{equation}
 where  $\mu$ is the dynamic molecular viscosity,
 $\kappa$ the volume viscosity (also called the second viscosity, usually neglected in the code, excepted for compressible flows).
$\tens{S}$ is called the strain rate and  $\Gamma$ is an eventual mass source term.

%-------------------------------------------------------------------------------
\subsection{Segregated solver: SIMPLEC}
A fractional step scheme is used to solve the mass and momentum equations
(see Chorin \cite{Chorin:1968}).

The first step (predictor step) provides predicted velocity
components: they are determined in a coupled way solving a $3 \ncell \times 3 \ncell $ system.
The mass equation is taken into account during the second step
(corrector step): a pressure Poisson equation is solved and the mass fluxes
at the cell faces are updated.

See the \doxygenfile{navstv_8f90.html}{programmers reference of the dedicated subroutine} for further details
on the SIMPLEC solver.

\subsubsection{Prediction step}
In this section, a predicted velocity field $ \widetilde{\vect{u}} $ is obtained by solving
the momentum equation of \eqref{eq:timstp:navier_stokes}

\begin{equation}\label{eq:timstp:prediction}
\begin{array}{rcl}
\displaystyle \rho \dfrac{ \vect{\widetilde{u}} - \vect{u}^{n} } { \Delta t}
+
\underbrace {
\gradt \,  \vect{\widetilde{u}}^{n+\theta} \cdot \left( \rho\vect{u} \right)
}_{
\text{convection}
}%end underbrace
&= & \underbrace{
\divv \left[ \left(\mu+\mu_t\right) \left( \gradt \, \vect{\widetilde{u}}^{n+\theta}
          + \transpose{\left(\gradt \, \vect{\widetilde{u}}^{n+\theta}\right)}
          - \frac{2}{3} \trace \left( \gradt \, \vect{\widetilde{u}}^{n+\theta}  \right)\right)\right]
}_{
\text{viscous term}
}% end underbrace
\\
 &-& \grad \left[ \left(P^*\right)^n \right]
 - \divv \left( \rho \deviator{\tens{R}}  + 2 \mu_T \deviator{\tens{S}}\right)
 + \left( \rho -\rho_0 \right) \vect{g} \\
&\displaystyle +&
\underbrace{
\Gamma \left( \vect{u}^{in}-\vect{\widetilde{u}} \right)
}_{
\text{Mass source term}
}%end underbrace
-
\underbrace{
\rho \tens{K} \, \vect{\widetilde{u}}
}_{
\text{Head loss}
}%end underbrace
 +
\underbrace{
\vect{ST}_{u}^{exp}+ \tens{ST}_{u}^{imp}\, \vect{\widetilde{u}}
}_{
\text{user source terms}
}%end underbrace
\end{array}
\end{equation}

For more details, see \S~ \ref{ap:navstv} and \ref{ap:predvv}. The \doxygenfile{predvv_8f90.html}{programmers reference
of the dedicated subroutine} can also be helpful.

\subsubsection{Correction step}
The predicted velocity has \emph{a priori} non-zero divergence. The second step
corrects the pressure by imposing the nullity\footnote{or the density time-variation
is the corresponding algorithm is chosen.}
of the stationary constraint for
the velocity computed at time instant ${t^{n+1}}$.
We then solve:
\begin{equation}\label{eq:timstp:c0_correction_step}
\left\{\begin{array}{ r c l}
\displaystyle \dfrac{(\rho\vect{u})^{n+1}-(\rho\vect{\widetilde{u}})^{n+1}}{\Delta t} &=&
-\grad{\delta P^{n+\theta}} ,\\
\displaystyle
\dive(\rho\vect{u})^{n+1} = \Gamma ,\\
\end{array}\right.
\end{equation}
% FIXME add head loss to the source term
where the pressure increment $\delta P^{n+\theta}$ is defined as:
\begin{equation}
\delta P^{n+\theta}=P^{n+\theta}-P^{n-1+\theta} .
\end{equation}

\begin{remark}
The  $\rho$ and $\mu$ quantities remain constant over the course of both
steps. If there has been any variation in the interval, their values will be
modified at the start of the next time step, after the scalars
(temperature, mass fraction,...) have been updated.
\end{remark}


For more details, see \S~ \ref{ap:navstv} and \ref{ap:resopv}. The \doxygenfile{resopv_8f90.html}{programmers reference
 of the dedicated subroutine} can also be helpful.

%-------------------------------------------------------------------------------
\subsection{Variable density conservative solvers}

Some flows such as those encountered in fire are natively unsteady.
In addition, density varies with mixture and temperature. Thus, the temporal
variation of density in Navier-Stokes equation has to be considered to get
unsteady phenomena such as flame puffings.
The SIMPLEC algorithm detailed before rewrites for a variable density field:

\begin{itemize}
 \item the prediction step solves the momentum equation with an explicit pressure:
%
\begin{equation}\label{eq:timstp:prediction}
\dfrac{\rho^n \vect{\widetilde{u}} - \rho^{n-1} \vect{u}^n}{\Delta t}
+ \divv \left[\vect{\widetilde{u}} \otimes (\rho \vect{u})^n \right] = -\grad p^n + \divv{\tens{\widetilde{\tau}}} + \rho^n\vect{g}.
\end{equation}

 \item the correction step computes the pressure increment $\delta p^{n+1} = p^{n+1} - p^n$
 used to correct the velocity respecting mass balance equation:
\begin{equation}
\left\{
  \begin{array}{r c l}
\dfrac{\rho^n \vect{u}^{n+1} - \rho^n\vect{\widetilde{u}}}{\Delta t} & = & -\grad \delta p^{n+1}\\
\dfrac{\rho^n - \rho^{n-1}}{\Delta t}+\divs{ \left(\rho \vect{u} \right)}^{n+1} & = & \Gamma
  \end{array}
\right.
\label{eq:timstp:correction}
\end{equation}
Combining the two equations of \eqref{eq:timstp:correction} leads to a Poisson equation on the pressure increment:
\begin{equation}
\divs{ \left( -\Delta t \grad \delta p^{n+1} \right) } = \Gamma
- \dfrac{\rho^n - \rho^{n-1}}{\Delta t} - \divs{ \left( \rho \vect{\widetilde{u}} \right) }
\end{equation}

\end{itemize}

%\begin{remark}
%Note that with this algorithm, the quantity $ $ is conserved
%\end{remark}

%-------------------------------------------------------------------------------
\subsection{Variable density semi-analytical solvers for fire modelling}
The algorithm detailed in this section is not conservative, and
is based on the state law for density.
The idea is to use an analytical expression for the total derivative of density,
expressed in terms of the solved variables
involved in the state law and the total derivatives of these scalars,
which are derived
from their non conservative balance equations ($\rho \DP{\varia} = \divs \left(K \grad Y \right)+ S_\varia $):
\begin{equation}
\dfrac{1}{\rho}\DP{\rho} = \dfrac{1}{\rho} \sum_{p=1}^N \der{\rho}{\varia_p}
\DP{\varia_p} = \dfrac{1}{\rho^2} \sum_{p=1}^N \der{\rho}{\varia_p} \left[\divs{ \left( K \grad {\varia_p} \right) }+S_{\varia_p}\right].
\label{eq:timstp:total_derivative}
\end{equation}

These scalars have to be independent to avoid counting several times
their effect on  the mass accumulation.
One has also to bear in mind that the total derivatives are dependent of the state law for density.

Therefore, the correction system \eqref{eq:timstp:correction} written in a non conservative form reads:
\begin{equation}
\left\{
  \begin{array}{rcl}
  \dfrac{(\rho \vect{u})^{n+1} - (\rho\vect{\widetilde{u}})}{\Delta t} = - \grad \delta p^{n+1} \\
  \DP{\rho} + \divs{ \left( \rho \vect{u} \right)^{n+1} } - \vect{u}^{n+1} \cdot \grad{\rho}^n = \Gamma,
  \end{array}
  \right.
\label{eq:timstp:correction3}
\end{equation}
and leads to a new equation for the pressure increment:
\begin{equation}
\DP{\rho} + \divs{ \left( \rho \vect{\widetilde{u}} \right) }
+ \divs{ \left( -\Delta t \grad \delta p^{n+1} \right) } - \vect{u}^{n+1} \cdot \grad{\rho}^n = \Gamma.
\label{eq:timstp:pression3}
\end{equation}
Using the expression for the corrected velocity:
\begin{equation}
\vect{u}^{n+1} = \vect{\widetilde{u}} - \dfrac{\Delta t}{\rho^n} \grad{\delta p^{n+1}},
\end{equation}
in equation \eqref{eq:timstp:pression3} leads to:
%
\begin{equation}
\dfrac{\Delta t}{\rho^n} \grad{\rho}^n \cdot \grad{\delta p^{n+1}}
- \divs{ \left( \Delta t \grad \delta p^{n+1} \right) } = \Gamma
+ \vect{\widetilde{u}} \cdot \grad{\rho}^n - \DP{\rho}
- \divs{ \left( \rho \vect{\widetilde{u}} \right) }.
\label{eq:timstp:pression4}
\end{equation}
%
Using some relationships between differential operators one obtains:
\begin{equation}
\dfrac{\Delta t}{\rho^n}\grad{\rho}^n \cdot \grad{\delta p^{n+1}}
= \divs{ \left( \dfrac{\Delta t}{\rho^n}\grad{\rho}^n \delta p^{n+1} \right) }
- \delta p^{n+1} \divs{ \left( \dfrac{\Delta t}{\rho^n} \grad{\rho^n} \right) },
\end{equation}
leads to a steady convection diffusion equation for the pressure increment:
\begin{equation}
\underbrace{ \divs{ \left( \dfrac{\Delta t}{\rho^n} \grad{\rho}^n \delta p^{n+1} \right) }
           - \divs{ \left( \dfrac{\Delta t}{\rho^n} \grad{\rho^n} \right) } \delta p^{n+1}
}_{\text{convection}}
-
\underbrace{
  \divs{ \left( \Delta t \grad \delta p^{n+1} \right) }
}_{\text{diffusion}}
 = \Gamma - \DP{\rho}
- \rho^n \divs{\vect{\widetilde{u}}}
\label{eq:timstp:pression5}
\end{equation}
This equation is solved in two steps with $\delta p = \delta p_1 + \delta p_2$
to split it into a pure diffusive step (which is the only remaining step when
the density is constant) and a convective/diffusive step:
%
\begin{enumerate}
 \item the first step solves the following equation which is closed to the equation solved in the standard algorithm,
\begin{equation}
 - \divs{ \left( \Delta t \grad \delta p_1^{n+1} \right) }
 = \Gamma - \DP{\rho} - \rho^n \divs{\vect{\widetilde{u}}}
 \label{eq:timstp:pression6}
\end{equation}
\item the second step solves the following steady convection diffusion equation
(see \S \ref{sec:algebrea_iterative_process}):
\begin{equation}
 \begin{array}{rcl}
  - \divs{ \left( \dfrac{\Delta t}{\rho^n} \grad{\rho^n} \right) } \delta p_2^{n+1}
  + \divs{ \left( \dfrac{\Delta t}{\rho^n} \grad{\rho}^n \delta p_2^{n+1} \right) }
  - \divs{ \left( \Delta t \grad \delta p_2^{n+1} \right) } = \\
    \divs{ \left( \dfrac{\Delta t}{\rho^n} \grad{\rho^n} \right) } \delta p_1^{n+1}
  - \divs{ \left( \dfrac{\Delta t}{\rho^n} \grad{\rho}^n \delta p_1^{n+1} \right) }
  \end{array}
 \label{eq:timstp:pression7}
\end{equation}
\end{enumerate}

%-------------------------------------------------------------------------------
%\subsection{Pseudo-coupled solver}

%\section{Steady algorithm}

%\section{Unsteady algorithm}