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

% This file is part of Code_Saturne, a general-purpose CFD tool.
%
% Copyright (C) 1998-2019 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.

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

\programme{predvv}\label{ap:predvv}

\hypertarget{predvv}{}

\vspace{1cm}
%-------------------------------------------------------------------------------
\section*{Function}
%-------------------------------------------------------------------------------
This subroutine implements the prediction step for the velocity field $\vect{u}$.
This consists of solving the momentum equation by treating the pressure explicitly.
Once the pressure correction step has been performed in the subroutine \fort{resopv},
the velocity-pressure solution is obtained using the mass conservation law:
\begin{equation}
\frac{\partial \rho } {\partial t}+ \dive(\rho \vect{u}) = \Gamma,
\end{equation}
where $\Gamma$ is the mass source term\footnote{In $kg.m^{-3}.s^{-1}$ }.\\
The Reynolds-averaged momentum conservation equation obtained by applying
the fundamental theorem of dynamics is:
\begin{equation}
\frac {\partial (\rho \vect {u})} {\partial t }+
\dive(\rho \vect{u} \otimes \vect{u}) =
\dive(\underline{\underline{\sigma}}) + \vect{S} - \dive{(\rho\,\tens{R})}\end{equation}
where:
\begin{equation}
\underline{\underline{\sigma}} = - p \underline{\underline{Id}} + \underline{\underline{\tau }}
\end{equation}
with the following linear relationship for the Newtonian flows:
\begin{equation}
\begin{array}{lcl}
&\displaystyle \underline{\underline{\tau}} = 2\ \mu\ \underline{\underline{D}}
+\,
 \lambda\ tr(\underline{\underline{D }})\ \underline{\underline{Id}} &\\
&\displaystyle \underline{\underline{D}}=\frac{1}{2}\ (\ggrad \underline
{u} +\ ^t\ggrad \vect {u})
\end{array}
\end{equation}


$\tens{\sigma}$ represents the stress tensor, $\tens{\tau}$ the viscous stress tensor,
$\mu$ the dynamic viscosity (molecular and potentially turbulent),$\tens{D}$
the rate of deformation tensor\footnote{Not to be confused, despite the same
notation $D$, with the diffusive fluxes described in the subroutine
\fort{navstv}}, $\tens{R}$ the Reynolds which appears when applying the
average operator to the simultaneous equation, $\vect{S}$ the source
terms.\\
$\lambda$ is the second coefficient of viscosity. It is related to the volume viscosity
 $\kappa$ via the expression
\begin{equation}
\lambda=\kappa-\frac{2}{3}\mu
\end{equation}
When the Stokes hypothesis holds, the volume viscosity $\kappa$ is
zero, in other words $3\lambda+2\mu=0$. This hypothesis is implicit in the code
and in the rest of this, except in the compressible module.\\


The equation for the conservation of momentum is finally written as:
\begin{equation}
\begin{array}{lcl}
&\displaystyle \rho\,
\frac{\partial \vect {u} } {\partial t} = -\
\underbrace {\dive(\rho \vect{u} \otimes \vect{u})}_{\text{
convection}} +\ \underbrace {\dive (\mu\ \ggrad \vect {u})}_{\text{
diffusion}} &\\
&\displaystyle \underbrace { +\ \dive (\mu \,^t\ggrad \vect {u}) }_{\text{
transpose of the velocity gradient term}}
\underbrace { - \ \frac {2} {3}\ \grad (\mu\ \dive \vect {u})}_{\text{
secondary viscosity}}\ \ - \dive{(\rho \tens{R})}
 -\ \grad(p) + (\rho -\rho_0)\,\vect {g} +
\vect{u}\,\dive (\rho\,\vect {u})&\\
&\displaystyle +\underbrace {\Gamma
(\vect{u}_{\,i}-\vect{u})}_{\text{mass source term}}-
\underbrace {\rho\
\tens{K} \vect {u}}_{\text{head loss}} +
\underbrace { \vect{T}_{\,s}^{\,exp}+
T_{s}^{\,imp}\ \vect{u}}_{\text{other source terms}}
\label{Base_Predvv_eqqdm}

\end{array}
\end{equation}
with $p$ denoting the difference in pressure relative to the reference hydrostatic
pressure (the real hydrostatic pressure being calculated with the density
$\rho$ and not with $\rho_{\,0}$):
\begin{equation}
p=p^*-\rho_{\,0}\ \vect{g}\cdot\vect{X}
\end{equation}
(\vect{X} being the vector of components $x$, $y$ and $z$).\\
$\mu_t$, $\tens{K}$, $\vect{u}_{\,i}$ represent respectively
the turbulent dynamic viscosity, the tensor related to head losses and the value of the
variable associated with the mass source.\\
The divergence of the Reynolds stress tensor is expressed by:
\begin{equation}
-\dive{(\rho\,\tens{R})}=
\begin{cases}
\vect{0} & \text{in laminar}, \\
 -\displaystyle\frac {2} {3}\, \grad (\mu_t\ \dive \vect {u})+\dive (\mu_t\ (\ggrad \vect {u}+ \,^t\ggrad \vect {u}))-\frac {2}{3}\,\grad (\rho\, k) & \text{for turbulent}\\
 & \text{viscosity}, \\
 -\dive(\rho\,\tens{R})& \text{for second-order}\\
 & \text{models},\\
-\displaystyle\frac {2} {3}\, \grad (\mu_t\ \dive \vect {u})+\dive (\mu_t\ (\ggrad \vect {u}+ \,^t\ggrad \vect {u})) & \text{in  LES}\\
\end{cases}
\end{equation}
The mass source term involves terms corresponding to the local velocity field
$\vect {u}$ as well as a velocity field $\vect {u}_{\,i}$
associated with the mass injected (or extracted). When $\Gamma<0$, we remove mass from the system thereby
obtaining $\vect{u}_{\,i} = \vect{u}$. The mass term is zero (\emph{i.e.} $\Gamma
(\vect{u}_{\,i}-\vect{u})= \vect{0} $). When $\Gamma>0$, a non-zero
mass term if $\vect{u}_{\,i} \ne \vect{u}$.
All terms appearing in the conservation of momentum equation, apart from
the convection and diffusion terms, are computed in this subroutine and
transmitted to the
subroutine \fort{codits} which solves the full equation (convection-diffusion with source terms).

See the \doxygenfile{predvv_8f90.html}{programmers reference of the dedicated subroutine} for further details.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discretization}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The convective term in $\dive(\vect{u} \otimes \rho\,\vect{u})$
introduces a nonlinearity and coupling of the components of the velocity
$\vect{u}$ in the equation (\ref{Base_Predvv_eqqdm}). The three velocity
components are linearized and decoupled during discretization of the momentum
equation in this prediction step.\\
For instance, taking:
\begin{equation}
\vect{\widetilde{u}}= \vect{u}^n + \delta \vect{u}
\end{equation}
The exact contribution of the convective term to take into account in the prediction step would be:
\begin{equation}\label{Base_Predvv_Conv_exact}
\begin{array}{ll}
\dive(\vect{\widetilde{u}} \otimes \rho\,\vect{\widetilde{u}}) =
\dive(\vect{u}^{n} \otimes \rho\,\vect{u}^{n}) + \dive(\delta \vect{u} \otimes
\rho\,\vect{u}^{n}) +  \underbrace { \dive(\vect{u}^{n} \otimes
\rho\,\delta \vect{u})}_{\text {linear coupling term}} +  \underbrace { \dive(\delta \vect{u} \otimes
\rho\,\delta \vect{u})}_{\text {nonlinear coupling term}}\\
\end{array}
\end{equation}
The last two terms of the expression (\ref{Base_Predvv_Conv_exact}) are {\em a priori} neglected
so as to obtain a velocity system that is decoupled and thereby avoid the
inversion of a potentially very large matrix. These two terms can
nevertheless be taken into account in an approximate manner by explicit
extrapolation of the mass flux in $n+\theta_F$ (for the linear coupling term
arising from the convection of $\vect{u}^{n}$ by $\delta \vect{u}$) and by
applying a fixed-point subiteration over the subroutine \fort{navstv} (for
the nonlinear term) implemented by specifying  $\var{NTERUP}>1$.

The equation (\ref{Base_Predvv_eqqdm}) is discretized at the time level
$n+\theta$ using a $\theta$-scheme and the tensor of velocity head losses
decomposed into the sum of its diagonal $\tens{K}_{d}$ and extra diagonal $\tens{K}_{e}$ terms (such that
 $\tens{K}=\tens{K}_{d}+\tens{K}_{e}$).\\
$\bullet$ The pressure is assumed to be known at the point in time $n-1+\theta$
(pressure/velocity time lag) and the gradient evidently calculated at this instant.\\
$\bullet$ The secondary viscosity and the gradient transpose source terms,
those coming from the turbulence model\footnote{with the exception of $\dive (\mu_t\ (\ggrad
\vect {u}))$}, $\rho\,\tens{K}_{\,e}\ \vect{u}$, $(\rho -\rho_0)
\vect {g}$ as well as $\vect{T}_{s}^{\,exp}$ and
$\Gamma\,\vect{u}_{\,i}$ are evaluated explicitly at time $n$ or else
extrapolated according to the time scheme selected for the physical properties and the
source terms.\\
$\bullet$ The source terms $\vect{u}\,\,\dive (\rho\,\vect {u})$,
$\Gamma\,\,\vect{u}$, $T_{s}^{\,imp}\,\,\vect{u}$ and
$-\rho\,\tens{K}_{\,d}\,\,\vect{u}$ are implicit and calculated at
instant $n+\theta$.\\
$\bullet$ The diffusion term  $\dive (\mu_{\,tot}\,\ggrad \vect{u})$ is implicitized: the velocity is taken at the instant $n+\theta$ and the viscosity either explicitized or extrapolated.\\
$\bullet$ Lastly, the convection term in $\dive(\,\vect{u} \otimes
(\rho\vect{u})\,)$ is treated implicitly: the resolved velocity
component is
taken at $n+\theta$ and the mass flux determined explicitly or extrapolated
at
$n+\theta_F$.

For clarification, unless explicitly stated otherwise the physical properties
$\Phi$ ($\rho,\,\mu_{tot},\,...$) and the mass flux $(\rho\vect{u})$
are  to be evaluated respectively at time steps $n+\theta_\Phi$ and
$n+\theta_F$, with $\theta_\Phi$ and $\theta_F$ dependent upon the specific
time schemes selected for these quantities\footnote{cf. \fort{introd}}.

The time-discrete form of the momentum equation (\ref{Base_Predvv_eqqdm}) then reads:

\begin{equation}\label{Base_Predvv_eq_di1}
 \begin{array}{c}
\displaystyle \rho\,\ \frac{ \vect {\widetilde{u}}^{n+1} -\vect {u}^{n} }
{\Delta t} + \dive(\,\vect{\widetilde{u}}^{n+\theta} \otimes (\rho\vect{u})\,) -\dive
(\mu_{\,tot}\,\ggrad \vect{\widetilde{u}}^{n+\theta}) =
\\
\displaystyle
 - \grad p^{n-1+\theta} + \dive (\rho\,\vect {u})\,\vect{\widetilde{u}}^{n+\theta} +(\Gamma\,\vect{u}_{\,i})^{n+\theta_S}-\Gamma^n\,\,\vect{\widetilde{u}}^{n+\theta}
\\
\begin{array}{c}
\displaystyle
- \rho\,\tens{K}_{\,d}^{n}\,\,\vect{\widetilde{u}}^{n+\theta} - (\rho\,\tens{K}_{\,e}\ \vect{u})^{n+\theta_S} + (\vect{T}_{s}^{\,exp})^{\,n+\theta_S} + T_{s}^{\,imp}\,\,\vect{\widetilde{u}}^{n+\theta}
\\
\displaystyle
+[\dive (\mu_{\,tot}\,^t\ggrad \vect {u})]^{n+\theta_S}-\frac {2} {3}[\,\grad (\mu_{\,tot}\,\dive \vect {u})]^{n+\theta_S} + (\rho -\rho_0) \vect {g}
 - (\vect{turb})^{n+\theta_S}
\end{array}
\end{array}
\end{equation}
where, for simplification, the following definitions have been set:
\begin{equation}
\mu_{\,tot}=
\begin{cases}
\mu+\mu_t & \text{for turbulent viscosity or LES models}, \\
\mu & \text{for second order models or in laminar regimes}
\end{cases} \
\end{equation}
\\
and:
\begin{equation}
\vect{turb}^{n}=
\begin{cases}
\displaystyle\frac {2}{3}\grad (\rho^{n}\,k^{n}) & \text{for turbulent viscosity models}, \\
\dive(\rho^{n}\,\tens{R}^n) & \text{for second order models},\\
0 & \text{in laminar or LES}\\
\end{cases}
\end{equation}
By analogy with the $\theta$-scheme equation for a scalar variable, $\,
\vect {\widetilde{u}}^{n+\theta}$ is interpolated from the predicted velocity
$\vect {\widetilde{u}}^{n+1}$ in the following manner\footnote{if
$\theta=1/2$, or a time extrapolation is used, second-order is only obtained provided the time step $\Delta t$ is constant and spatially uniform.}~:
\begin{equation}
\vect {\widetilde{u}}^{n+\theta}=\theta\, \vect
{\widetilde{u}}^{n+1}+(1-\theta)\, \vect {u}^{n}\\
\end{equation}
With:
\begin{equation}
\left\{
\begin{array}{ll}
\theta = 1   & \text{For a first-order implicit Euler scheme.}\\
\theta = 1/2 & \text{For a second-order Crank-Nicolson scheme.}\\
\end{array}
\right.
\end{equation}

Using the above interpolation function, (\ref{Base_Predvv_eq_di1}) is then rewritten in the following form:

\begin{equation}\label{Base_Predvv_eq_di2}
\begin{array}{c}
\displaystyle \underbrace{\left(\frac{\rho}{\Delta t} -\theta \,\dive (\rho\,\vect {u}) +\theta \,\, \Gamma^n +
\theta \,\, \rho\,\tens{K}_{\,d}^n-\theta \,T_s^{\,imp} \right)}_{\displaystyle f_s^{imp}}\, (\vect {\,\widetilde{u}}^{n+1} -\vect {u}^{n})
\\
 +\, \theta\, \dive(\vect {\widetilde{u}}^{n+1} \otimes (\rho\vect{u}))-\, \theta\,\dive (\mu_{\,tot}\,\ggrad \vect {\widetilde{u}}^{n+1}) =
\\
-\,(1-\theta)\, \dive(\vect {u}^{n} \otimes (\rho\vect{u})) +\,(1-\theta)\,\dive (\mu_{\,tot}\,\ggrad \vect {u}^{n})
\\
f_s^{exp}\left\{
\begin{array}{c}
\displaystyle
- \grad p^{n-1+\theta} + \dive (\rho\,\vect {u})\,\vect{u}^{n} +\,(\,\Gamma^{n}\,\vect{u}_{\,i}\,)^{n+\theta_S}- \Gamma^n\,\,\vect{u}^{n}
\\
\displaystyle
-(\,\rho\,\tens{K}_{\,e}\ \vect{u}\,)^{n+\theta_S} -\rho\,\tens{K}_{\,d}^n\ \vect{u}^{n}+ (\vect{T}_{s}^{\,exp})^{\,n+\theta_S} + T_s^{\,imp}\,\,\vect {u}^{n}
\\
\displaystyle
+[\dive (\mu_{\,tot}\,^t\ggrad \vect {u}\,)]^{n+\theta_S}-\frac {2} {3}[\,\grad (\mu_{\,tot}\,\dive \vect {u}\,)]^{n+\theta_S} + (\rho -\rho_0) \vect {g}-(\vect{turb})^{n+\theta_S}
\end{array}
\right.
\end{array}
\end{equation}

from which the equation solved by the subroutine \fort{codits} is obtained:
\begin{equation}\begin{array}{c}
\displaystyle
f_s^{\,imp}(\vect {\widetilde{u}}^{n+1}-\vect {u}^{n}) + \theta\, \dive(\vect{\widetilde{u}}^{n+1} \otimes (\rho
\vect{u})) - \theta\,\dive (\,\mu_{\,tot}\,\ggrad \vect{\widetilde{u}}^{n+1}) =
\\\\
\displaystyle
-(1-\theta)\,\dive(\vect{u}^{n} \otimes (\rho \vect{u}))+(1-\theta)\,\dive (\,\mu_{\,tot}\,\ggrad \vect{u}^{n})
+ \vect{f}_{\,s}^{\,exp}
\end{array}
\end{equation}
The spatial discretization scheme is elaborated further in the section relating to the subroutine \fort{codits}.\\



\minititre{Remarks:}
{\tiny$\blacksquare$} In the standard case without extrapolation, the term
$-\,T_s^{\,imp}$ is only added to the coefficient $f_s^{\,imp}$ if it is positive so as not to weaken
the diagonal dominance of the matrix to invert.\\
{\tiny$\blacksquare$} If, on the other hand, an extrapolation is used,
$\,T_s^{\,imp}$ is added to $f_s^{\,imp}$ whatever its sign. In effect, the intuitive idea which consists of taking:
\begin{equation}
\begin{cases}
\displaystyle
(\vect{T}_{s}^{\,exp} + T_{s}^{\,imp}\,\vect {u})^{\,n+\theta_S} &
\text{if } T_{s}^{\,imp} > 0\\
\displaystyle
(\vect{T}_{s}^{\,exp})^{\,n+\theta_S} + T_{s}^{\,imp}\,\vect{u}^{n+\theta} &\text{if not}\\
\end{cases}
\end{equation}
leads to inconsistency (loss of coherence) in the numerical treatment of the problem should $T_s^{imp}$ changes sign between two time steps.\\
{\tiny$\blacksquare$} The diagonal part $\tens{K}_{\,d}$ of the head loss tensor term is used in $f_s^{\,imp}$. Strictly speaking, only the positive contributions should in fact be retained (as pointed out in the associated user subroutine \fort{uskpdc}). The  approach currently in use requires improvement.\\
{\tiny$\blacksquare$} The term $\theta\,\Gamma^{n}-\theta\,\dive
(\rho\,\vect {u})$ is not problematic for the diagonal dominance of the matrix as it is exactly counterbalanced by the convection term (cf. \fort{covofi}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The conservation of momentum equation is solved in a decoupled manner. The integration of the different terms has therefore been carried out in order to treat separately the equations obtained for each component of the velocity.\\
For each of these components, the second member $f_s^{exp}$ of the system  (\ref{Base_Predvv_eq_di2}), the implicit terms of the system (with the exception of the convection-diffusion terms) and the term representing the total viscosity at the  internal\footnote{value required for the integration of the diffusion term in \fort{codits}, $\displaystyle(\mu_{\,tot})_{ij}\frac{\var{SURFN}}{\var{DIST}}$} and boundary faces is calculated in the subroutine \fort{predvv}. These terms are then transmitted to the subroutine \fort{codits} which constructs and then solves the resulting complete system of equations, with the convection-diffusion terms, for each component of the velocity.
\\\\
The normalized residual for convergence of the solution of the system for the pressure correction (\fort{resopv}) is calculated in
\fort{predvv}. It is defined as the norm of magnitude
  $$\dive(\rho\,\widetilde{\vect{u}}^{n+1}+\Delta t \grad{P^{n-1+\theta}})-\Gamma$$
integrated over each cell \var{IEL} of the mesh ($\Omega_{iel}$) or, symbolically, as the square root of the sum over the mesh cells of the quantity
\begin{center}
\var{XNORMP(IEL)}=
$\int\limits_{\Omega_{iel}}[{\dive(\rho\,\widetilde{\vect{u}}^{n+1}+\Delta\,t\,\grad{P^{n-1+\theta}})
-\Gamma\,]\,d\Omega}$.
\end{center}

It represents the second member of the equation system that would be solved for the pressure if the pressure gradient were not taken into account in the velocity prediction step. Note that, if the second member of the pressure increment equation were to be used directly, a normalized residual tending to zero would be obtained for a stationary calculation that had been run to convergence. This result would be penalizing and of little use for the computations.

At the start of \fort{predvv}, $\widetilde{\vect{u}}^{n+1}$ is not yet available and it is not possible therefore to calculate the total normalized residual. On the other hand, calculating the total residual at the end of \fort{predvv} is undesirable as this would require monopolising a working array to store the pressure gradient for the duration of the subroutine \fort{predvv}.
The calculation of the normalized residual is therefore carried out in two stages.

The quantity $\int\limits_{\Omega_{iel}}{\dive(\Delta\,t\,\grad{P^{n-1+\theta}}) -\Gamma d\Omega}$ is calculated at the beginning of \fort{predvv} and the complement $\int\limits_{\Omega_{iel}}{\dive(\rho\,\widetilde{\vect{u}}^{n+1})d\Omega}$ added at the end of \fort{predvv}.

The pressure gradient at the cells is thus first computed at instant $n-1+\theta$ by a call to \fort{grdcel}. The subroutine \fort{inimas} is then used to evaluate $\Delta\,t\,S\,\grad{P^{n-1+\theta}}\cdot\vect{n}$ at the faces (of surface $S$ and normal $\vect{n}$). On entry into \fort{inimas}, the working array \var{TRAV} contains $\frac{\Delta\,t}{\rho}\,\grad{P^{n-1+\theta}}$; when exiting this subroutine, the \var{VISCF} and \var{VISCB} arrays contain the value of $\Delta\,t\,S\,\grad{P^{n-1+\theta}}\cdot\vect{n}$ at the internal and boundary faces respectively.

We then use \fort{divmas} takes the value of $\int\limits_{\Omega_{iel}}{\dive(\Delta\,t\,\grad{P^{n-1+\theta}}) d\Omega}$ at the cells from the contents of \var{VISCF} and \var{VISCB} and places it into the array \var{XNORMP}.

Lastly, the contribution $\int\limits_{\Omega_{iel}}{ -\Gamma^{n} d\Omega}$ of the mass source term is added to \var{XNORMP}.

For $\rho\,\widetilde{\vect{u}} + \Delta t \grad{P}$, we apply the velocity boundary conditions. For the computation of
$\int\limits_{\Omega_{iel}}{\dive(\Delta\,t\,\grad{P^{n-1+\theta}}) d\Omega}$, the boundary conditions for the pressure gradient (or more precisely for $\frac{\Delta\,t}{\rho}\,\grad{P^{n-1+\theta}}$) are therefore the homogenized boundary conditions of the velocity. The pressure gradient (or rather $\frac{\Delta\,t}{\rho}\,\grad{P^{n-1+\theta}}$) is consequently assumed to be zero normal to the walls and the inlets and to remain constant in the direction normal to symmetries and to the outlets.

In addition, to save computation time when passing through \fort{inimas}, the face values on non-orthogonal meshes are evaluated merely to first-order accuracy in space (no reconstruction: \var{NSWRP=1}). Local accuracy is moreover of no particular interest as the aim here is simply to determine a normalized global residual.

The computation of the residual will be completed at the end of \fort{predvv}.


\etape{Partial calculation of the normalized residual for the pressure step}

During this first step, we use the array \var{XNORMP(NCELET)} to compute the quantity $$\dive(\Delta t \,\grad{P^{n-1+\theta}})-\Gamma$$ integrated over each cell \var{IEL} of the mesh ($\Omega_{iel}$), which is expressed as:
$$\var{XNORMP(IEL)}=\int\limits_{\Omega_{iel}}{\dive(\Delta\,t\,\grad{P^{n-1+\theta}})-\Gamma
d\Omega}$$
This operation is carried out by successively implementing \fort{inimas}
(calculation at the faces in \var{VISCF} and \var{VISCB} of $\Delta t\,\grad{P^{n-1+\theta}}$ using the contents of the working array \var{TRAV}=$\frac{\Delta t}{\rho} \grad{P^{n-1+\theta}}$ together with homogeneous velocity boundary conditions and without reconstruction) followed by \fort{divmas} (computation in \var{XNORMP} of the integral over the cells).  With a simple loop, we then add the mass source term contribution $\Gamma$.
This calculation is completed at the end of \fort{predvv}.\\

\etape{Partial calculation of the term $\vect{f}_s^{\,exp}$}

To represent the second member corresponding to each component of the velocity, we use the arrays \var{TRAV}(\var{IEL},\var{DIR}), \var{TRAVA}(\var{IEL},\var{DIR}) and \var{PROPCE}, with \var{IEL} being the number of the cell and \var{DIR} the direction (x, y, z). Four cases need to be considered depending on whether or not the source terms are extrapolated at $n+\theta_S$ and whether or not a fixed-point iteration is implemented on the velocity-pressure system (\var{NTERUP}$>1$).
\\
$\bullet$ If the source terms are extrapolated and we iterate over \fort{navstv}\\
\begin{itemize}
\item [-]\var{TRAV} receives the source terms which are recomputed during each iteration over \fort{navstv} but are not extrapolated
($\grad{P^{n-1+\theta}}$ and $(\rho -\rho_0) \vect {g}$\footnote{In reality, the value of
$(\rho -\rho_0) \vect {g}$ doesn't change, however this is a computation and it avoids the need for additional treatment of this term.}).\\
\item [-]\var{TRAVA} receives the source terms whose values do not change during the course of the iterations over \fort{navstv} nor are they extrapolated
($T_s^{imp}\,u^n\,,-\rho\,\tens{K}_{d}\,u^n\,,\,-\Gamma^n\,u^n,...$).\\
\item [-]\var{PROPCE} receives the source terms that have to be extrapolated.\\
\\
\end{itemize}
$\bullet$ When there is no iteration over \fort{navstv}, \var{TRAVA} has no use and its contents are directly stored in \var{TRAV}.\\
\\
$\bullet$ If there is no source term extrapolation, \var{PROPCE} serves no purpose and its contents are directly stored in \var{TRAVA} (or in \var{TRAV} if \var{TRAVA} is also of no use).\\
Therefore, when there is neither source term extrapolation nor iteration over \fort{navstv}, all the source terms go directly into \var{TRAV}.\\
\\
\begin{itemize}
\item We already have a pressure gradient on the cells at the instant in time $n-1+\theta$. The gravity term is then added to the vector \var{TRAV} which already contains pressure gradient. This means we have, for example, for the direction x:
\begin{equation}
\var{TRAV}(\var{IEL},1) = |\Omega_{IEL}| (-\displaystyle (\frac {\partial p}
{\partial x})_{\var{IEL}}+(\rho(\var{IEL})-\rho_0)g_x)
\end{equation}

\item If source term extrapolation is used, at the first iteration over \fort{navstv}, the vector \var{TRAV} (or \var{TRAVA}) receives $-\theta_S$ times the contribution at time $n-1$ of the source terms to be extrapolated\footnote{Because
$(\vect{T}_s^{exp})^{n+\theta_S}=(1+\theta_S)\,(\vect{T}_s^{exp})^n
-\,\theta_S\, (\vect{T}_s^{exp})^{n-1}$} (stored in \var{PROPCE}). \var{PROPCE} is then reinitialised to zero so as to be able to receive afterwards the contribution at the current time step of the extrapolated source terms.

\item The value of the term corresponding to the turbulence model is only calculated during the first iteration over \fort{navstv} then added to \var{TRAVA}, \var{TRAV} or \var{PROPCE} depending on whether the source terms are extrapolated, and/or whether we iterate over \fort{navstv}.\\

{\tiny$\blacksquare$} Turbulent viscosity models:\\
If $\var{IGRHOK}=1$, then we calculate $-\displaystyle \frac{2}{3}\ \rho\ \grad k$ (and not, as should be the case,
$-\displaystyle \frac{2}{3}\grad (\rho k)$) by way of simplification
(cf. paragraph~\ref{Base_Predvv_section4}). The gradient of the turbulent kinetic energy $k$ is computed on the cell by the subroutine \fort{grdcel}. \\
If $\var{IGRHOK}=0$, this term is expected to be implicitly taken into account in the pressure.\\
{\tiny$\blacksquare$} Second order models:\\
Computation of the term $-\dive(\rho \tens{R})$ is implemented in two steps. A call is first made to the subroutine \fort{divrij}, which projects the vector $\tens{R}.\vect{e}_{\var{DIR}}$ onto the cell faces along the direction \var{DIR}, following which we then call the subroutine \fort{divmas} to compute the divergence.\\
\linebreak
\item The secondary viscosity $- \displaystyle \frac {2} {3}\grad (\mu_{\,tot} \dive \vect\,{u})$ and transposed gradient $ \dive (\mu_{\,tot} \,^t\ggrad \vect {u})$ terms, when these taken into account (\emph{i.e.} \var{IVISSE}\,(\var{IPHAS})\ = 1, where \var{IPHAS} is the number of the phase treated) are computed by the subroutine \fort{visecv}. They are only calculated at the first iteration over \fort{navstv}. During this step, the array \var{TRAV} serves as a working array when the subroutine  \fort{visecv} is being called. It reverts to its normal value at the end of the call during which its contents are stored temporarily in the vectors \var{W7} to \var{W9}.
\\
\item The terms corresponding to the velocity head losses ($\rho \tens{K}
{u}$), if there are any (\ $\var{NCEPDP}> 0$\ ), are calculated by the subroutine \fort{tspdcv} at the first iteration over \fort{navstv}. They are decomposed into two parts:\\
{\tiny$\blacksquare$} The first part, corresponding to the contribution of the diagonal terms ($-\rho\,\tens{K}_{\,d}\vect{u}$), is not extrapolated.\\
{\tiny$\blacksquare$} The second, which corresponds to the extra diagonal terms
($-\rho\,\tens{K}_{\,e}\vect{u}$), may or may not be extrapolated.\\
During this step, the array \var{TRAV} is used as a working array while the call is being made to the subroutine \fort{tspdcv}. It reverts to its normal value at the end of this call, its contents having been meanwhile stored in the vectors \var{W7} to \var{W9}.

\end{itemize}

\etape{Calculation of the viscosity term at the faces
$\displaystyle(\mu_{\,tot})_{ij}\frac{\var{SURFN}}{\var{DIST}}$}
The calculation of the total viscosity term at the faces is performed by the subroutine \fort{viscfa} and the solutions stored in the arrays \var{VISCF} and \var{VISCB} for the internal and boundary faces respectively.\\
During integration of the convection-diffusion terms in the subroutine \fort{codits}, the non-reconstructed terms, integrated in the matrix $\tens{EM}$, are singled out from the full set of terms (unreconstructed +
reconstruction gradients) associated to the operator $\mathcal{E}$ (nonlinear)\footnote{ Ensuring coherence with the definition of the operators $\mathcal{EM}$ and $\mathcal{E}$ employed in \fort{navstv} }.
Likewise, we differentiate between the total viscosity at the faces used in $\mathcal{E}$, arrays \var{VISCF} and \var{VISCB}, and the total viscosity at the faces employed in $\tens{EM}$, arrays \var{VISCFI} and \var{VISCBI}.\\
For turbulent viscosity models and in LES, these two arrays are identical and contain $\mu_t+\mu$.
For second order models, they normally contain $\mu$ however, for the simple reason of numerical stability, we can choose to put $\mu_t+\mu$ in the matrix (\textit{i.e.} in $\tens{EM}$) while keeping $\mu$ in the second member (\textit{i.e.} in $\mathcal{E}$). Because the solution is obtained by increments, this manipulation doesn't change the result. This option is activated by the indicator $\var{IRIJNU}\ =\ 1$\\
When there is no diffusion of the velocity (\ \var{IDIFF}(\var{IUIPH})\ $<$\ 1), the terms \var{VISCF} and \var{VISCB} are set to zero.
\linebreak

\etape{Calculation of the whole second member, of $f_s^{\,imp}$ and solving of the equation}
The evolution equations for the components of the momentum are solved in a coupled manner. Consequently, we use a single array, \var{ROVSDT}, to represent the diagonal of the matrix obtained for each of the velocity components.\\
For each of these components:\\

\begin{itemize}
\item During the first iteration over the subroutine \fort{navstv}, the implicit and explicit parts of the user-prescribed source terms are computed by a call to the subroutine \fort{ustsns}.\\
{\tiny$\blacksquare$} The implicit part ($T_s^{imp}$) is saved in the vector \var{XIMPA} for subsequent iterations if the fixed-point method is used on the velocity-pressure system, with the contribution resulting from the same implicit terms ($T_s^{imp}\,\vect{u}^n$) being added either to \var{TRAVA} or to \var{TRAV}.\\
{\tiny$\blacksquare$} The explicit part ($T_s^{exp}$) is added to \var{TRAVA}, \var{TRAV} or \var{PROPCE} depending on whether the source terms are extrapolated or not, and whether we iterate over \fort{navstv}.\\

\item The mass accumulation term ($\dive(\rho \vect {u})$) is calculated by calling the subroutine \fort{divmas} with the mass flux as argument. During the first iteration of the subroutine \fort{navstv},
the term corresponding to the explicit contribution of the mass accumulation ($\vect{u}^{n}\ \dive(\rho \vect{u})$) is added to \var{TRAVA} or to \var{TRAV}. The vector \var{ROVSDT} is initialized by $\theta\,\dive(\rho \vect{u})$ (to remain coherent with the approach implemented in the subroutine
\fort{bilsc2}) after which the contribution of the non-stationary term ($\displaystyle
\frac{\rho}{\Delta t}$) is added to the latter.\\

\item The vector \var{ROVSDT} is next completed with the contribution of the user-defined implicit source terms (stored in \var{XIMPA}) and with that of the head loss term ($\rho\,\tens{K}_{\,d}$) should $\var{NCEPDP}>0$.\\
{\tiny$\blacksquare$} When there is no source term extrapolation, the implicit part of the user source terms is  only added to \var{ROVSDT} if it is negative so as not to weaken the diagonal of the system.\\
{\tiny$\blacksquare$} When the source terms are extrapolated however, it is taken into account whatever its sign may be.\\

\item The implicit and explicit mass source terms, if any (~$\var{NCESMP}>0$~), are computed at the first iteration by the subroutine  \fort{catsma} over \fort{navstv}. $\Gamma\,\vect{u}_i$ is added to \var{TRAV}, \var{TRAVA} or \var{PROPCE} for potential extrapolation. $\Gamma\,\vect{u}^n$ is added to \var{TRAV} or \var{TRAVA} and
$-\Gamma$ to \var{ROVSDT}.
\\
\item The second member is then assembled taking account of all of the
contributions stored in the arrays \var{PROPCE}, \var{TRAVA} and
\var{TRAV}.\\
{\tiny$\blacksquare$} If the source terms are extrapolated, then the second member reads:
$$\var{SMBR}=(1-\theta_S)\,\var{PROPCE}+\var{TRAVA}+\var{TRAV}$$
{\tiny$\blacksquare$} If not, we directly obtain:
$$\var{SMBR}=\var{TRAVA}+\var{TRAV}$$

\item If a specific physics (lagrangian, electric arc, ...) is taken into account, its contribution is directly added to \var{SMBR}.
\\
\item The solution of the linear system is computed by the subroutine
\fort{codits}, taking as argument \var{ROVSDT} and \var{SMBR}.\\

\item When strengthened transient velocity-pressure coupling ($\ \var{IPUCOU} = 1\ $) (only available with a first-order scheme, where there is neither source term extrapolation nor iteration over \fort{navstv}) is implemented, \fort{codits} solves the following equation:
\begin{equation}\label{Base_Predvv_Eq_Tdir}
\tens{EM}_{\,\var{DIR}}\,.\, (\tens{RHO}^{\,n})^{-1}\,.\,\vect{T}_{\,\var{DIR}} =
\tens{\Omega}\,.\,\vect{1}
\end{equation}
with $\tens{RHO}^n$ the diagonal tensor of element $\rho^{n}_{IEL}$,
$\tens{\Omega}$ the diagonal tensor of element $|\Omega_{IEL}|$ and $\vect{1}$ denoting the
vector whose components are all equal to one.\\
Inversion of the system by \fort{codits} yields
$(\tens{RHO}^{\,n})^{-1}\,.\,\vect{T}_{\,\var{DIR}}$, which is subsequently multiplied by $\tens{RHO}^{\,n}$
to obtain $\vect{T}_{\,\var{DIR}}$.
This process is repeated for each component \var{DIR} of the velocity. $\vect{T}_{\,\var{DIR}}$
is therefore a diagonal matrix approximation of
$\tens{RHO}^{\,n}\,.\,\tens{EM}_{\,\var{DIR}}^{-1}$, with
$\tens{EM}_{\,\var{DIR}}$ representing the implicit part of the solution to the discrete momentum equation (\emph{i.e.} \var{ROVSDT} + contribution of the convection-diffusion terms accounted for in the subroutine
\fort{matrix}). $\vect{T}_{\,\var{DIR}}$ is also involved in the correction step (cf. subroutine \fort{resopv}).\\
\end{itemize}
This terminates the loop over the velocity components.\\

\etape{Final calculation of the normalized residual for the pressure step}

As mentioned beforehand, the computation of the normalized residual for the pressure step of \fort{resopv} can now be completed.

The \var{XNORMP} array already contains
$\int\limits_{\Omega_{iel}}{\dive(\Delta\,t\,\grad{P^{n-1+\theta}}) -\Gamma d\Omega}$ to which we now add
$\int\limits_{\Omega_{iel}}{\dive(\rho\,\widetilde{\vect{u^{n+1}}})d\Omega}$.

To do so, we follow the same procedure used previously to compute the integral
$\int\limits_{\Omega_{iel}}{\dive(\Delta\,t\,\grad{P^{n-1+\theta}}) d\Omega}$. A call to
\fort{inimas} enables to obtain
$\rho\,S\,\widetilde{\vect{u}}^{n+1}\cdot\widetilde{\vect{n}}$ at the faces calculated from the of
$\widetilde{\vect{u}}^{n+1}$ known at the cells (array \var{RTP}). The boundary conditions
employed in \fort{inimas} are naturally those of the velocity. As beforehand, evaluation of the face values is only up to first-order in space on non-orthogonal meshes (no reconstruction~: \var{NSWRP=1}) in order to save computation time when passing through \fort{inimas}. The subroutine \fort{divmas} is then used to compute the divergence $\int\limits_{\Omega_{iel}}{\dive(\rho\,\widetilde{\vect{u}}^{n+1})d\Omega}$ at the cells, which is then added directly to \var{XNORMP}.

Lastly, the normalized residual is determined and then stored in \var{RNORMP(IIPHAS)} by a call to \fort{prodsc}  (which computes the sum of the squares of the values at the faces contained in \var{XNORMP} and then takes the square root of the result).\\


\newpage

The different contributions (excluding convection-diffusion) assigned to each of the vectors \var{TRAV}, \var{TRAVA}, \var{PROPCE} and \var{ROVSDT} at iteration $n$ are summarized in tables (\ref{Base_Predvv_tab_ext0}), (\ref{Base_Predvv_tab_ext1}), (\ref{Base_Predvv_tab_exp0}) and (\ref{Base_Predvv_tab_exp1}). We differentiate two cases for each of the time-stepping schemes applied to the source terms, depending on whether or not a fixed-point algorithm is used to solve the velocity-pressure system (iteration over \fort{navstv} for
\var{NTERUP}$>1$). Unless otherwise specified, the physical properties $\Phi$
($\rho$,$\mu$,etc...) are assumed to be evaluated at the time level $n+\theta_\Phi$, and the mass flux
 $(\,\rho \vect{u})$ at $n+\theta_F$, where $\theta_\Phi$ and $\theta_F$ are dependent on the specific time-stepping schemes selected for these quantities (cf. \fort{introd}).
\\\\
The terms appearing in these tables are written as they have been implemented in the code, this being the source of some differences in relation to the form adopted in equation (\ref{Base_Predvv_eq_di2}).
\\\\
In order to simplify the problem, we pose:
\begin{equation}\notag
\mu_{\,tot}=
\begin{cases}
\mu+\mu_t & \text{for turbulent viscosity models and in LES}, \\
\mu & \text{for second-order models and in laminar regime}\\
\end{cases} \
\end{equation}

\minititre{With extrapolation of source terms}
\begin{equation}\notag
\vect{turb}^{n}=
\begin{cases}
\displaystyle\frac {2}{3}\,\rho^{n}\,\grad (\,k^{n}) &
\text{for turbulent viscosity models}, \\
\dive(\rho^{n}\,\tens{R}^n) & \text{for second-order models},\\
0 & \text{in laminar or LES computations.}\\
\end{cases}
\end{equation}
$\bullet$ \var{NTERUP} $=$ 1 : $\var{SMBR}^n=(1-\theta_S)\,\var{PROPCE}^n+\var{TRAV}^n$
\begin{equation}\label{Base_Predvv_tab_ext0}
\begin{array}{|l|c|}
\hline
\var{ROVSDT}^{n}
& \displaystyle
\frac{\rho}{\Delta t} -\theta \,\dive (\rho\,\vect {u}) +\theta \,\, \Gamma^n + \theta \,\, \rho\,\tens{K}_{\,d}^n-\theta \,T_s^{\,imp} \\
\hline
\var{PROPCE}^{n}
& \displaystyle
\vect{T}_{s}^{\,exp,\,n}-\,\rho^{n}\,\tens{K}_{\,e}^{n}\ \vect{u}^{n} + \,\Gamma^{n}\,\vect{u}_{\,i}^{n}\\
& \displaystyle
-\vect{turb}^{n}+ \dive (\mu_{\,tot}^{n}\,^t\ggrad \vect {u}^{n}\,)+ \frac {2} {3}\,\grad (\mu_{\,tot}^{n}\,\dive \frac{(\rho \vect {u})}{\rho^{n}}\,)\\
\hline
\var{TRAV}^{n} & \displaystyle
- \grad p^{n-1+\theta} + (\rho -\rho_0) \vect {g} \\
& \displaystyle
-\theta_S\,\var{PROPCE}^{n-1} -\rho\,\tens{K}_{\,d}^n\ \vect{u}^{n} \\
& \displaystyle
+ T_s^{\,imp}\,\,\vect {u}^{n} + \dive (\rho\,\vect {u})\,\vect{u}^{n} - \Gamma^n\,\,\vect{u}^{n}\\
\hline
\end{array}
\end{equation}
\\\\
$\bullet$ \var{NTERUP} $>$ 1 (sub-iteration $k$) : $\var{SMBR}^n=(1-\theta_S)\,\var{PROPCE}^n+\var{TRAVA}^n+\var{TRAV}^n$
\begin{equation}\label{Base_Predvv_tab_ext1}
\begin{array}{|l|c|}
\hline
\var{ROVSDT}^{n}
& \displaystyle
\frac{\rho}{\Delta t} -\theta \,\dive (\rho\,\vect {u}) +\theta \,\, \Gamma^n + \theta \,\, \rho\,\tens{K}_{\,d}^n-\theta \,T_s^{\,imp} \\
\hline
\var{PROPCE}^{n}
& \displaystyle
\vect{T}_{s}^{\,exp,\,n}-\,\rho^{n}\,\tens{K}_{\,e}^{n}\ \vect{u}^{n} + \,\Gamma^{n}\,\vect{u}_{\,i}^{n}\\
& \displaystyle
-\vect{turb}^{n}+ \dive (\mu_{\,tot}^{n}\,^t\ggrad \vect {u}^{n}\,)+ \frac {2} {3}\,\grad (\mu_{\,tot}^{n}\,\dive \frac{(\rho \vect {u})}{\rho^{n}}\,)\\
\hline
\var{TRAVA}^{n} &
\displaystyle
-\theta_S\,\var{PROPCE}^{n-1} -\rho\,\tens{K}_{\,d}^n\ \vect{u}^{n} + T_s^{\,imp}\,\,\vect {u}^{n} + \dive (\rho\,\vect {u})\,\vect{u}^{n} - \Gamma^n\,\,\vect{u}^{n}\\
\hline
\var{TRAV}^{n}
& \displaystyle
- \grad (p^{n+\theta})^{(k-1)} + (\rho -\rho_0) \vect {g} \\
\hline
\end{array}
\end{equation}

\minititre{Without source term extrapolation}

\begin{equation}\notag
\vect{turb}^{n}=
\begin{cases}
\displaystyle\frac {2}{3}\,\rho\,\grad (\,k^{n}) &
\text{for turbulent viscosity models}, \\
\dive(\rho\,\tens{R}^n) & \text{for second-order models},\\
0 & \text{in laminaire or LES computations.}\\
\end{cases}
\end{equation}
$\bullet$ \var{NTERUP} $=$ 1 : $\var{SMBR}^n=\var{TRAV}^n$
\\\\
\begin{equation}\label{Base_Predvv_tab_exp0}
\begin{array}{|l|c|}
\hline
\var{ROVSDT}^{n} &
\displaystyle \frac{\rho}{\Delta t} -\theta \,\dive (\rho\,\vect {u}) +\, \Gamma^n + \, \rho\,\tens{K}_{\,d}^n+Max(-\,T_s^{\,imp},0)\\
\hline
\var{TRAV}^{n}
& \displaystyle
- \grad p^{n-1+\theta} + (\rho -\rho_0) \vect {g} \\
& \displaystyle
+ \vect{T}_{s}^{\,exp}-\,\rho\,\tens{K}_{\,e}^{n}\ \vect{u}^{n} + \,\Gamma^{n}\,\vect{u}_{\,i}^{n}\\
& \displaystyle
-\vect{turb}^{n}+ \dive (\mu_{\,tot}\,^t\ggrad \vect {u}^{n}\,)+ \frac {2} {3}\,\grad (\mu_{\,tot}\,\dive \frac{(\rho \vect {u})}{\rho}\,)\\
& \displaystyle
-\rho\,\tens{K}_{\,d}^n\ \vect{u}^{n}+ T_s^{\,imp}\,\,\vect {u}^{n} + \dive (\rho\,\vect {u})\,\vect{u}^{n} - \,\Gamma^{n}\,\vect{u}^{n}\\
\hline
\end{array}
\end{equation}
\\\\
$\bullet$ \var{NTERUP} $>$ 1 (sub-iteration $k$) : $\var{SMBR}^n=\var{TRAVA}^n+\var{TRAV}^n$
\begin{equation}\label{Base_Predvv_tab_exp1}
\begin{array}{|l|c|}
\hline
\var{ROVSDT}^{n} &
\displaystyle \frac{\rho}{\Delta t} -\theta \,\dive (\rho\,\vect {u}) +\, \Gamma^n + \, \rho\,\tens{K}_{\,d}^n+Max(-\,T_s^{\,imp},0)\\
\hline
\var{TRAVA}^{n} &
\displaystyle
\vect{T}_{s}^{\,exp}-\,\rho\,\tens{K}_{\,e}^{n}\ \vect{u}^{n} + \,\Gamma^{n}\,\vect{u}_{\,i}^{n}\\
& \displaystyle
-\vect{turb}^{n}+ \dive (\mu_{\,tot}\,^t\ggrad \vect {u}^{n}\,)+ \frac {2} {3}\,\grad (\mu_{\,tot}\,\dive \frac{(\rho \vect {u})}{\rho}\,)\\
& \displaystyle
-\rho\,\tens{K}_{\,d}^n\ \vect{u}^{n}+ T_s^{\,imp}\,\,\vect {u}^{n} + \dive (\rho\,\vect {u})\,\vect{u}^{n} - \Gamma^n\,\,\vect{u}^{n}\\
\hline
\var{TRAV}^{n} &
\displaystyle
- \grad (p^{n+\theta})^{(k-1)} + (\rho -\rho_0) \vect {g} \\
\hline
\end{array}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points to treat}\label{Base_Predvv_section4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\etape{Evaluation of the term $\grad(\rho k)$ for turbulent viscosity models}
When modelling turbulent viscosity, we compute $\rho\ \grad k$ instead of $\grad(\rho k)$. This
approximation was implemented from the outset because the $\rho k$ boundary conditions are not directly accessible, in contrast to those of $k$.\\

\etape{Expression of $\tens{EM}$}
With incremental solving, it is not absolutely necessary for convergence that the value of the viscosity appearing in the expression of the operator $\mathcal{E}$ be the same as that taken into account in the incremental system matrix, $\tens{EM}$. In $R_{ij}-\varepsilon$, for example, the total viscosity used in $\tens{EM}$ contains the molecular viscosity but can also contain the turbulent viscosity if the option \var{IRIJNU = 1 } is selected, whereas only the former is included in the operator $\mathcal{E}$. This addition of the turbulent viscosity, which is not at all relevant to the $R_{ij}-\varepsilon$ model, has been inherited from practices implemented in
ESTET and N3S-EF to improve numerical stability (incremental smoothing). However, it may have other, potentially less desirable effects. Moreover, it has not been demonstrated to date that this practice is of absolute necessity in the use of \CS. An in-depth study would therefore be of some interest.


\etape{Normalized residual relating to the pressure step}
We could check the normalized residual computation and in particular the use of the velocity boundary conditions.

\etape{Calculation of the head losses}
When the source terms are extrapolated in time, we now obtain:
\begin{equation}\notag
(\tens{K}_{\,e}\vect{u})^{n+\theta_S} + \tens{K}_{\,d}^{n}\ \vect{u}^{n+\theta}
\end{equation}
It would also be possible to envisage using:
\begin{equation}\notag
(\tens{K}_{\,e}\vect{u})^{n+\theta_S} + \tens{K}_{\,d}^{n+\theta_S}\ \vect{u}^{n+\theta}
\end{equation}