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

% 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.

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

\programme{resopv}\label{ap:resopv}

\hypertarget{resopv}{}

\vspace{1cm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Function}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The velocity projection (or pressure correction) step is effected in this subroutine, called from \fort{navstv}. The equation of motion (prediction) is solved in \fort{predvv} with a fully explicit treatment of the pressure term. The resulting velocity field does not satisfy the continuity equation. Two correction algorithms are proposed:
\begin{enumerate}
\item The algorithm that we will call "weak velocity-pressure coupling". This algorithm is implemented extensively in industrial source codes. It only couples the velocity and the pressure through the mass term (it is the algorithm proposed by default). It consists of a \textit{SIMPLEC}-type algorithm, similar to \textit{SIMPLE} although the latter  accounts for the simplified diagonals of the convection, diffusion and implicit source terms in addition to the mass term.
\item The strengthened velocity-pressure coupling algorithm (option \var{IPUCOU  = 1}). This algorithm couples the velocity and pressure through all of the terms (convection, diffusion and implicit source terms) of the equation of motion, though it is still approximate. In practice, the advantage of this algorithm is that it allows large time steps without entirely decoupling the velocity and the pressure.
\end{enumerate}

Taking $\delta p$ as the pressure increment ({\it i.e.} $p^{n+1} = p^n+\delta p$) and $\widetilde{\vect{u}}$ the velocity field resulting from the prediction step, from a continuous point of view the projection step essentially comes down to solving a Poisson equation for the pressure:

\begin{equation}
  \dive(\,{\tens{T}^n\ \grad{\delta p}}) = \dive(\,{\rho \,\widetilde{\vect{u}}})
\end{equation}

and then correcting the velocity:

\begin{equation}
\vect{u}^{n+1} = \vect{u}^n - \frac{1}{\rho}\ \tens{T}^n\ \grad{\delta p}
\end{equation}

$\tens{T}^n$ is a second-order tensor whose components are homogeneous over a time step.

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discretization}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The velocity prediction step is described in \fort{predvv}. An operator notation is adopted here to provide a simplified explanation of the code-based algorithms.
The discrete equation of motion at an intermediate point in time of solution is written, in each direction of space $\alpha$ ($\alpha \in \{1,2,3\}$)in the form:

\begin{equation}
\displaystyle
\tens{M}_{\,\alpha}^n\ \tens{R}^n\ (\widetilde{\vect{V_\alpha}} - \vect{V_\alpha}^n) +
\tens{A_\alpha}^n\ \widetilde{\vect{V_\alpha}} =
-\ \tens{G_\alpha}\ \vect{P}^n + \vect{S_\alpha}^n + \tens{I}_{\,s,\alpha}\ \widetilde{\vect{V_\alpha}}
\label{Base_Resopv_QDM}
\end{equation}
\begin{itemize}
\item[$\star$] $\tens{M}_{\,\alpha}^n$, is a diagonal matrix, of dimension $ \var{NCEL} \times \var{NCEL}$, containing the ratio of the volume of a cell to the local time step ($\ \displaystyle {\tens{M}_{\,\alpha}^n}(i,i) = \frac{|\Omega_i|}{{\Delta t}_{\,\alpha,I}^n} $\ ), ${\Delta t}_{\alpha,I}^n$ represents the time step at time level $n$ in the (spatial) direction $\alpha$ at the cell $\Omega_i$,
\item[$\star$] $ \tens{R}^n $, is the diagonal matrix of dimension $\var{NCEL} \times \var{NCEL}$, that contains the density (which is separated out from the mass matrix during this step so that a projection can be made of the value of $\rho \vect{u} $). By definition, $ \tens{R}^n(i,i) = \rho^{\,n}_{\,I}$, which means that the appearance of a vacuum (or null space) is naturally excluded and the matrix will therefore always be invertible,
\item[$\star$] $\widetilde{\vect{V_\alpha}}$, of dimension $\var{NCEL}$, is an array in which the $\alpha^{\text{th}}$ component of the predicted velocity field $\widetilde{\vect{u}}$ is stored,
\item[$\star$] $\vect{V_\alpha}^n$, of dimension $\var{NCEL}$, is an array storing the $\alpha^{\text{th}}$ component of the velocity $\vect{u}^n$ at the previous instant of time $n$,
\item[$\star$] $\tens{A_\alpha}^n$ denotes the convection/diffusion operator (it is not necessarily linear owing to the possible use of slope tests in the convection scheme and may be dependent on $\vect{V_\alpha}$),

\item[$\star$] $\tens{G_\alpha}$ is the linear, "cell" gradient operator\footnote{strictly speaking, this operator would no longer be truly linear if a gradient constraint option were to have been activated by the user} in the direction $\alpha$ (it is therefore applied on the vectors of dimension \var{NCEL}),
\item[$\star$] $\vect{P}^n$, of dimension $\var{NCEL}$, is an array used to store the pressure $p^n$ computed at each cell during the previous time step,

\item[$\star$] $\vect{S_\alpha}^{\,n}$ is the array of dimension $\var{NCEL}$ that contains all of the explicit source terms (see \fort{predvv} for more detail),

\item[$\star$] $\tens{I}_{\,s,\alpha}$ is the diagonal tensor related to the implicit source terms of the velocity components.
\end{itemize}

The correction step consists of imposing the continuity constraint for the mass conservation:
\begin{equation}\label{Base_Resopv_Eq_Cont}
\dive(\,{\rho\,\vect{u}}\,) = \Gamma
\end{equation}
where $\Gamma$ is an eventual mass source term.
%($\Gamma = 0$ for incompressible and dilatable flows).

Let $\vect W$ be the array of dimension $3 \times \var{NCEL}$ that contains all of the momentum components ($\vect{V}$ denotes $\vect{V}^n$, $\vect{V}^{n+1}$ or $\widetilde{\vect{V}}$).
$$ \vect{W} = \tens{R}^n\ \vect{V} = \left(
                    \begin{array} {c}
                    \rho^n\ \vect{V_1} \\
                    \rho^n\ \vect{V_2} \\
                    \rho^n\ \vect{V_3}
                    \end{array}
              \right)
$$

Let $\tens{D}$ be the divergence operator.
The continuity equation (\ref{Base_Resopv_Eq_Cont}) is rewritten in compact form as:
\begin{equation}\notag
\tens{D}\ \vect{W} = \vect{\Gamma}
\end{equation}
$\vect{\Gamma}$ contains the values of $\Gamma$ at the cell centroids.

Rearranging the discrete equation (\ref{Base_Resopv_QDM}), we write for all $\alpha \in \{1,2,3\}$ :
\begin{equation}
\displaystyle
(\tens{M}_{\,\alpha}^n + \tens{A_{\,\alpha}}^n\ {(\tens{R}^{n})}^{-1} -\ \tens{I}_{\,s,\alpha}\ {(\tens{R}^n)}^{-1})\ \tens{R}^n\ \widetilde{\vect{V_\alpha}} = -\ \tens{G_\alpha}\ \vect{P}^n + \vect{S_{\,\alpha}}^n + \tens{M_{\,\alpha}}\ \tens{R}^n\ \vect{V_\alpha}^n
\label{Base_Resopv_eqn2}
\end{equation}

By grouping and posing the following definitions:
$$ \tens{B_{\,\alpha}} =  \tens{M}_{\,\alpha}^n + \tens{A_{\,\alpha}}^n\ {(\tens{R}^n)}^{-1} - \tens{I}_{\,s,\alpha}\ {(\tens{R}^n)}^{-1}$$
$$ \tens{B} = \left(
                    \begin{array} {ccc}
                    \tens{B_1} & 0 & 0 \\
                    0 & \tens{B_2} & 0 \\
                    0 & 0 & \tens{B_3}
                    \end{array}
              \right)
$$
$$ \tens{G} = \left(
                    \begin{array} {c}
                    \tens{G_1} \\
                    \tens{G_2} \\
                    \tens{G_3}
                    \end{array}
              \right)
$$
$$ \vect{S}^{\,n} = \left(
                      \begin{array} {c}
                      \vect{S_1}^n+\ \tens{M}_{\,1}\ \tens{R}^n\ \vect{V_1}^n\\
                      \vect{S_2}^n+\ \tens{M}_{\,2}\ \tens{R}^n\ \vect{V_2}^n\\
                      \vect{S_3}^n+\ \tens{M}_{\,3}\ \tens{R}^n\ \vect{V_3}^n
                      \end{array}
                \right)
$$

we can thus write the simplified equation as:
\begin{equation}
\tens{B}\ \widetilde{\vect{W}} = -\ \tens{G}\ \vect{P}^{\,n} + \vect{S}^{\,n}
\label{Base_Resopv_eqn3}
\end{equation}

With the fractional step method, (\ref{eq:continuousnavstv}) is decomposed into a sequence of two steps:
\begin{enumerate}
\item solution of the equation (\ref{Base_Resopv_QDM}) (this equation yields equation(\ref{Base_Resopv_eqn3})), namely:
\begin{equation}
\displaystyle
\tens{M}_{\,\alpha}^n\ \tens{R}^n\ (\widetilde{\vect{V_\alpha}} - \vect{V_\alpha}^n) +
\tens{A_{\,\alpha}}^n\ \widetilde{\vect{V_\alpha}} - \tens{I}_{\,s,\alpha}\ \widetilde{\vect{V_\alpha}}
 =
-\ \tens{G_\alpha}\ \vect{P}^{\,n} + \vect{S_\alpha}^n
\label{Base_Resopv_QDM2}
\end{equation}

\item subtraction\footnote{We neglect any variation in the explicit source terms $\vect{S_\alpha}^n$ which are still those estimated at $n$.} of the solved velocity prediction equation (\ref{Base_Resopv_QDM2}) from the equation of motion evaluated at the subsequent time level $(n+1)$:
\begin{equation}
\displaystyle
\tens{M}_{\,\alpha}^n\ \tens{R}^n\ (\vect{V_\alpha}^{n+1} - \widetilde{\vect{V_\alpha}}) +
\tens{A_{\,\alpha}}^n\ (\vect{V_\alpha}^{n+1} - \widetilde{\vect{V_\alpha}}) - \tens{I}_{\,s,\alpha}\ (\vect{V_\alpha}^{n+1} - \widetilde{\vect{V_\alpha}})
 =
-\ \tens{G_{\,\alpha}}\ (\vect{P}^{\,n+1}- \vect{P}^{\,n})
\label{Base_Resopv_frac2}
\end{equation}
\end{enumerate}

Reverting to compact notation, Equation (\ref{Base_Resopv_frac2}) gives:

\begin{equation}
\tens{B}\ (\vect{W}^{n+1}-\widetilde{\vect{W}}) = -\ \tens{G}\ (\vect{P}^{\,n+1}-\vect{P}^{\,n})
\label{Base_Resopv_eqn5}
\end{equation}

with:
$$ \vect{W}^{\,n+1} = \tens{R}^n\ \vect{V}^{\,n+1}
$$

The continuity constraint has yet to be enforced:
$$\tens{D}\ \vect{W}^{\,n+1} = \vect{\Gamma} $$

By combining the equations of continuity and motion and postulating that $ \delta \vect{P} = \vect{P}^{\,n+1} - \vect{P}^{\,n}$, the following Poisson-type equation is derived:

\begin{equation}
 \tens{D}\ \tens{B}^{-1}\,\tens{G}\ \,\delta \vect{P} =\ \tens{D}\ \vect{\widetilde{W}} -\ \vect{\Gamma}
\label{Base_Resopv_eqn5bis}
\end{equation}

We still need to invert the system (\ref{Base_Resopv_eqn5bis}) in order to determine $\delta p$ (and thus $p^{\,n+1}$) and correct the projected velocity field so as to obtain $\vect{u}^{\,n+1}$. The velocity correction is handled in \fort{navstv}, by incrementing the velocity by the calculated magnitude of the gradient of the pressure increment $\delta p$.

The problem initially arises with the calculation of $\tens{B}^{-1}$. It has already been judged too expensive computationally to calculate the inverse of $\tens{B}$. The solutions computed by the "weak velocity-pressure coupling" and the "strengthened velocity-pressure coupling" algorithms correspond to an approximation of this operator.

In the case of the "weak velocity-pressure coupling" algorithm, we assume $\tens{B}_{\,\alpha}^{-1} = \tens{M}_{\,\alpha}^{-1}$ (we could also include the diagonal terms of the convection, diffusion and implicit source terms).

With "strengthened velocity-pressure coupling", we invert the system\footnote{$\vect{\Omega} = (|\Omega_1|,...,|\Omega_{\var{NCEL}}|,|\Omega_1|,...,|\Omega_{\var{NCEL}}|,|\Omega_1|,...,|\Omega_{\var{NCEL}}|)$.} $\tens{B}\ \vect{T} = \vect{\Omega} $ and we define the equality $\tens{B}_{\,\alpha}^{-1}~=~\text{diag}(T_\alpha) $. This step takes place in the subroutine \fort{predvv}.

The use of operators when writing out the equations presents a major inconvenience when used in conjunction with collocated discretization. More specifically, the operator\footnote{We emphasize again that the operator $\tens{G}$ is the "cell gradients" operator applied to the explicit pressure during velocity prediction.}
$\tens{D}\ \tens{B}^{-1}\ \tens{G}$ can lead to odd-even decoupling of the nodes on a regular Cartesian
mesh\footnote{If $u_i$ is the value of a variable at the cell centres on a 1D Cartesian mesh, the Laplacian
of $u$ calculated by this operator in $i$ is written: $\displaystyle
\frac{u_{i-2}+2u_i-u_{i+2}}{4h^2}$, where $\displaystyle h$ is the space step.  This is where the decoupling of the cells originates.}.
To avoid this problem, we use the operator $\tens{L}$ (already present in the collocated finite volume
formulation of the operator $\dive(\,\grad\,))$ defined in each cell\footnote{Recalling that $Neigh(i)$ is the set of cell centres of the neighbouring cells of ${\Omega_i}$ and
$\gamma_b(i)$ the set of centres of the boundary faces, if any, of ${\Omega_i}$.} $\Omega_i$ of centre $I$ by\footnote{If $u_i$ is the value of a variable at the cell centres on a 1D Cartesian mesh, the Laplacian
of $u$ calculated by the latter operator in $i$ reads: $\displaystyle \frac{u_{i-1}+2u_i-u_{i+1}}{4h^2}$, with $\displaystyle h$ the spatial step.}:
\begin{equation}\label{Base_Resopv_eqn7}
\begin{array}{ll}
&(\tens{L}\ \delta \vect{P})_I = \sum\limits_{j\in Neigh(i)}[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}
+ \sum\limits_{k\in {\gamma_b(i)}}[\ \tens{T}^{\,n}_{\,b_{\,ik}}\ (\grad{\delta p})_{\,f_{\,b_{\,ik}}}]\,.\,\vect{S}_{\,b_{\,ik}}
\end{array}
\end{equation}
$\tens{T}^n$ is a diagonal tensor of order 2 containing the time steps in the three spatial directions with $\vect{S}_{\,ij}$ and $\vect{S}_{\,b_{\,ik}}$ the surface vector respectively of the purely internal face and of the boundary face $ik$. The gradient $(\grad{ \delta p})_{\,f_{  }}$  present in equation (\ref{Base_Resopv_eqn7}) is a facet\footnote{On orthogonal mesh, $\displaystyle (\grad{p})_{\,f_{\,ij}} \ . \ \vect{S}_{\,ij} = \frac{ p_J - p_I }{\overline{IJ}} S_{\,ij}$. $I$ (resp. $J$) et $I'$ (resp. $J'$) are in effect superimposed.} gradient.

From a continuous perspective, we can write\footnote{In the "weak velocity-pressure coupling" algorithm case, $\widetilde{T}^n_I = \Delta t^{\,n}_I$.}:
$$(\tens{L}\ \vect{P}^{n+1})_I = \int_{\Omega_i}\dive{(\ {\tens{T}}^n \ \grad{p^{n+1}})} \ d\Omega $$
%$$(\tens{L}\ \vect{P}^{n+1})_I = \int_{\Omega_i}\dive{(\ \widetilde{\tens{T}}^n \ \grad%{p^{n+1}})} \ d\Omega $$
The operator $\tens{L}$ does not pose a problem of odd/even decoupling on regular Cartesian meshes.

With the "weak velocity-pressure coupling" algorithm:
$ \tens{T}_{\,I}^{\,n} = \left(
                      \begin{array} {ccc}
                      \Delta t_I^{\,n} & 0 & 0 \\
                      0 & \Delta t_I^{\,n} & 0 \\
                      0 & 0 & \Delta t_I^{\,n}
                      \end{array}
                \right) $\\
and with the "strengthened velocity-pressure coupling" algorithm:
$ \tens{T}_{\,I}^{\,n} = \left(
                      \begin{array} {ccc}
                      T_{\,11,\,I}^{\,n} & 0 & 0 \\
                      0 &  T_{\,22,\,I}^{\,n}&0 \\
                      0 & 0 & T_{\,33,\,I}^{\,n}
                      \end{array}
                \right) $

The pressure system matrix is not easy to compute when "strengthened velocity-pressure coupling" is being used, this being due to the term:
$$[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta
p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}$$ for an internal face $ij$ and the term: $$[\ \tens{T}^{\,n}_{\,b_{\,ik}}\ (\grad{\delta p})_{\,f_{\,b_{\,ik}}}]\,.
\,\vect{S}_{\,b_{\,ik}}$$ for a boundary face.\\
The difficulty lies with the viscosity, which changes according to the spatial direction. As the tensor  $\tens{T}^n$ is effectively anisotropic, the pressure gradients for each direction must consequently be computed as a function of the normal gradient.
%(\delta p) = \tens{T}^{\,n}\ (\grad{(\delta p)})_{\,f } =

Without reference, for the moment, to the nature of the gradient $\grad$, we define for every scalar $a$:
$$ \vect{\widetilde{G}}\ a = \tens{T}^{\,n}\ \grad{a} =
              \left(
                    \begin{array} {c}
                    \displaystyle T_{11}^{\,n}\ \frac{\partial a}{\partial x}\\
                    \displaystyle T_{22}^{\,n}\ \frac{\partial a}{\partial y}\\
                    \displaystyle T_{33}^{\,n}\ \frac{\partial a}{\partial z}
                    \end{array}
              \right)
$$
We can also define:\\
$$ \left[(\vect{\widetilde{G}}\ a )_{\,cell}\right]_{\,ij}\,.\,\vect{S}_{\,ij} \overset{\text{\it\small def}}{=} \left[\tens{T}^{\,n}\ (\grad{a})\right]_{\,ij}\,.\,\vect{S}_{\,ij} $$  where
$\grad{a}$ is the standard cell gradient of $a$.\\

and likewise for the facet gradient $(\grad{a})_{\,f}$ of $a$:
$$ \left[(\vect{\widetilde{G}}\ a )_{\,f}\right]_{\,ij}\,.\,\vect{S}_{\,ij} \overset{\text{\it\small def}}{=} \left[\tens{T}^{\,n}\ ((\grad{a}))_{\,f}\right]_{\,ij}\,.\,\vect{S}_{\,ij} $$
\\

We need to calculate $\vect{\widetilde{G}}\ (\delta p)\,.\,\vect{S}$ at the face.\\
This is relatively simple with a cell gradient as all the components of the latter are perfectly calculable.\\
The difficulty lies, on the contrary, with a facet gradient because the only exploitable options is the decomposition of the face normal gradient\footnote{According to the formula $ \displaystyle\frac { p_{\,J'} - p_{\,I'}}{\overline{I'J'}} S_{\,ij}$.}. We need the quantities $\displaystyle \frac{\partial (\delta p)}{\partial x}$, $\displaystyle \frac{\partial (\delta p)}{\partial y}$ and $\displaystyle \frac{\partial (\delta p)}{\partial z}$ both explicitly and separately, which makes it difficult to use the normal pressure increment gradient.\\
We therefore make an approximation of the gradient of the pressure increment by assuming it to be equal to its normal component\footnote{{\it i.e.} no account is taken of the tangential component $((\grad{(\delta p)})_{\,f}\,.\,\vect{\tau}).\vect{\tau}$, $\vect{\tau}$ denoting the tangential unit vector.}, which is:
\begin{equation}
(\grad{(\delta p)})_{\,f} \approx ((\grad{(\delta p)})_{\,f}\,.\,\vect{n})\ \vect{n}
\end{equation}
where $\vect{n}$ denotes the outer unit normal vector.
We then obtain:
$$ \vect{\widetilde{G}}\ (\delta p) \approx ((\grad{(\delta p)})_{\,f}\,.\,\vect{n})\ (\tens{T}^n\ \vect{n})$$
This enables reducing the computations to those of a scalar time step $\widetilde{T}^n_{\,ij}$, given by:
\begin{equation}
\widetilde{T}^n_{\,ij} = (\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{n}
\label{Base_Resopv_approx2}
\end{equation}
In which case:
\begin{equation}
\begin{array} {lll}
[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}& \approx (\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{S}_{\,ij}\ (\grad{(\delta p)}_{f_{\,ij}}\,.\,\vect{n}) \\
 & =  (\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{n} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{n}) \ S_{\,ij} \\
 & =  \widetilde{T}^n_{\,ij} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{n}) \ S_{\,ij}
 \end{array}
\label{Base_Resopv_approximation}
\end{equation}
or on the other hand:
\begin{equation}
\begin{array} {lll}\label{Base_Resopv_Eq_Exacte}
\widetilde{T}^n_{\,ij} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{n}) \ S_{\,ij} &= \widetilde{T}^n_{\,ij} \ {\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{S}_{\,ij}\\
&=\,\widetilde{T}^n_{\,ij}\ \left[ P_J - P_I
 + (\ \vect{JJ'} - \vect{II'}\ )\,\displaystyle \frac{1}{2}\ (\,{\grad P}_I + \,{\grad P}_J\,)\right]\, \displaystyle \frac{S_{\,ij}}{\overline{I'J'}}
\end{array}
\end{equation}
This is used during reconstruction of the gradients in the second member of the final system of equations to solve, by calling the subroutines \fort{itrmas} and \fort{itrgrp}.\\
The approximation ($\displaystyle \widetilde{T}^n_{\,ij}=(\tens{T}^n_{\,ij}\ \vect{n})\,.\,\vect{n}$) is not in fact used. This is because the directional gradient of the pressure increment is computed using the subroutine \fort{grdcel}, allowing us to take the liberty of using the tensor $\tens{T}$ to correct the term $[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}$ .\\
In practice, the latter is discretized, for an internal face $ij$, by\footnote{The factor $\displaystyle \frac{1}{2}$ appearing in the discretization (\ref{Base_Resopv_grad1}) is introduced for reasons of numerical stability. Although the use of weighting coefficients on the faces would in effect yield a more accurate numerical solution, it would probably also be less stable.}
:
\begin{equation}\label{Base_Resopv_grad1}
\begin{array} {lll}
&[\ \tens{T}^{\,n}_{\,ij}\ (\grad{\delta p})_{\,f_{\,ij}}]\,.\,\vect{S}_{\,ij}
&=\left[\widetilde{T}^n_{\,ij}\ \displaystyle \frac{P_J - P_I}{\overline{I'J'}} + \displaystyle
 \frac{\vect{JJ'} -\vect{II'}}{\overline{I'J'}} \, \frac{1}{2}\ (\,\tens{T}^{\,n}_{\,I}\ {\grad P}_I + \tens{T}^{\,n}_{\,J}\ {\grad P}_J\,)\right]\,S_{\,ij}
\end{array}
\end{equation}
rather than by (\ref{Base_Resopv_Eq_Exacte}).
We use the same approach for the boundary terms.\\
The expression $\widetilde{T}^n_{\,ij} \ ({\grad{(\delta p)}}_{f_{\,ij}}\,.\,\vect{S}_{\,ij})$ will nevertheless be seen subsequently.\\
The last issue relating to the inversion of the system (\ref{Base_Resopv_eqn5bis}) concerns the tensor term  $\tens{D}\ \widetilde{\vect{W}}$. This term tends, \textit{via} the cell pressure gradients present in the equation of motion, to decouple the odd and even-numbered cells on a Cartesian grid.
To avoid this problem, we apply a variant of the Rhie \& Chow filter which enables to dissipate (or smooth) the pressure field contribution in the equation of motion.
Expressed in discretized form, this yields:
\begin{equation}\label{Base_Resopv_eqn8}
\begin{array}{lll}
(\tens{D}\ \widetilde{\vect{W}})_I = &\sum\limits_{j\in Neigh(i)}[\rho^{\,n} \widetilde{\vect{u}} + \alpha_{\,Arak}\ (\vect{\widetilde{G}}\ (p^n))_{\,cell}]_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} - \alpha_{\,Arak}\sum\limits_{j\in Neigh(i)}\ \widetilde{T}^n_{\,ij}\ (\grad{p^{n}})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} \\
&+ \sum\limits_{k\in {\gamma_b(i)}}[\rho^{\,n} \widetilde{\vect{u}} + \alpha_{\,Arak}\ (\vect{\widetilde{G}}\ (p^n))_{\,cell}]_{\,f_{\,b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}} - \alpha_{\,Arak}\sum\limits_{k\in {\gamma_b(i)}}\widetilde{T}^n_{\,b_{\,ik}}\ (\grad{p^{n}})_{\,f_{b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}}\\
&= \sum\limits_{j\in Neigh(i)} \widetilde{m}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{m}_{\,b_{\,ik}}
\end{array}
\end{equation}
where $\widetilde{m}_{\,ij}$ and $\widetilde{m}_{\,b_{\,ik}}$ denote respectively the mass flux at the purely internal face and the boundary face $ik$) modified by the Rhie \& Chow filter.\\
$(\grad{p^{n}})_{f_{\,ij}}$ represents a facet gradient while $[\rho^{\,n} \widetilde{\vect{u}} + \alpha_{\,Arak}\ (\vect{\widetilde{G}}\ (p^n))_{\,cell}]_{\,f_{\,ij}}$ denotes a face value interpolated from the estimated cell values (the pressure gradient in this term is a cell-based gradient). This clarification also applies to the boundary terms.\\
For historical reasons, $ \alpha_{\,Arak} $ is known as "the Arakawa coefficient d'Arakawa" in \CS. It is denoted \var{ARAK}.

It should be reiterated that, for the Rhie \& Chow filter, the tensor $\tens{T}^n$ is used in the volumetric term $\widetilde{\vect{G}}$ whereas the approximation (\ref{Base_Resopv_approx2}) is applied when computing the facet gradient.

In accordance with equations (\ref{Base_Resopv_eqn5bis}) and (\ref{Base_Resopv_eqn7}), we finally solve the Poisson equation under the form:
\begin{equation}
\sum\limits_{j\in Neigh(i)}\ \widetilde{T}^n_{\,ij}\ (\grad{(\delta p)})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{T}^n_{\,b_{\,ik}}\ (\grad{(\delta p)})_{\,f_{b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}} = \sum\limits_{j\in Neigh(i)} \widetilde{m}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{m}_{\,b_{\,ik}} -\ \Gamma_I
\label{Base_Resopv_Poisson}
\end{equation}
To take account of the non-orthogonal elements, we use an incremental method to solve the linear equation system, including the reconstructed terms in the second member. If $\delta(\delta p)$ denotes the increment of the increment (the increment of the variable $\delta p$ to compute) and $k$ the fixed-point iteration index, we solve exactly :
\begin{equation}
\begin{array}{lcl}
&\sum\limits_{j\in Neigh(i)}\ \widetilde{T}^n_{\,ij}\ \displaystyle \frac{(\delta(\delta p))_I^{\,k+1}-(\delta(\delta p))_J^{\,k+1}}{\overline{I'J'}}\ S_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{T}^n_{\,b_{\,ik}}\ \displaystyle\frac{(\delta(\delta p))_I^{\,k+1}-(\delta(\delta p))_{\,b_{ik}}^{\,k+1}}{\overline{I'F}}\ S_{\,b_{\,ik}} \\
&= \sum\limits_{j\in Neigh(i)} \widetilde{m}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}\widetilde{m}_{\,b_{\,ik}}\\
& - \sum\limits_{j\in Neigh(i)}\,\widetilde{T}^n_{\,ij}\ (\grad{(\delta{p})^k})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} - \sum\limits_{k\in {\gamma_b(i)}}\,\widetilde{T}^n_{\,b_{\,ik}}\ (\grad{(\delta{p})^k})_{\,f_{b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}}
- \ \Gamma_I\\
&= \sum\limits_{j\in Neigh(i)} m^{\,k}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}}m^{\,k}_{\,b_{\,ik}} - \ \Gamma_I
\end{array}
\label{Base_Resopv_equation resolue}
\end{equation}
with:
\begin{equation}\label{Base_Resopv_increment}
\left\{\begin{array}{lll}
(\delta(\delta p))^{0} &= 0\\
(\delta(\delta p))^{k+1} &= (\delta p)^{k+1}-(\delta p)^k &\text {$\forall k \in \mathbb{N}$}\\
\end{array}\right.
\end{equation}
The facet gradients, denoted $(\grad{(\delta{p})^k})_{\,f_{\,ij}}$ and $(\grad{(\delta{p})^k})_{\,f_{b_{\,ik}}}$, will be reconstructed.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
We present hereafter the algorithm as it is written in \fort{resopv}.

$\vect{T}^n$ designates an array of dimension $3$ containing the local time steps in each direction (for use with ''strengthened velocity-pressure coupling"). We keep the same notation for the ''weak velocity-pressure coupling" algorithm although in this case the time steps are equal in the three spatial directions.

\etape {Computation of the matrix of the equation system to be solved}
\begin{itemize}
\item calculation of the diffusion coefficient at the cell faces for use in the Laplacian of pressure (the diffusion coefficient uses the calculation time step or that of the ''strengthened velocity-pressure coupling" algorithm). One of two cases will apply, depending on the algorithm selected by the user for the velocity-pressure coupling:

\begin{enumerate}
\item Call to \fort{viscfa} with a total viscosity equal to the time step $\Delta t_I^n$ for the "weak velocity-pressure coupling" algorithm (\var{IPUCOU} = 0),
\item Call to \fort{visort} with a diagonal total viscosity for the "strengthened velocity-pressure coupling" algorithm (\var{IPUCOU} = 1). It is at this level that $\widetilde{T}^n_{\,ij}$ is calculated. The equivalent time steps calculated in the subroutine \fort{predvv} beforehand are stored in the array \var{TPUCOU}.
\end{enumerate}

\item Call to \fort{matrix} for construction of the pressure diffusion matrix (without the reconstruction terms which cannot be taken into account if we wish to preserve a sparse matrix structure) using the previously calculated viscosity and the array \var{COEFB} of pressure boundary conditions $p^n$ (we impose a homogeneous Neumann condition on $\delta p$ for a Neumann condition on $p$ and \textit{vice versa}.

\end{itemize}

\etape{Calculation of the normalized residual \var{RNORMP}}
Since version 1.1.0.s, this step is accomplished within the subroutine \fort{predvv} and the normalized residual transmitted \textit{via} the variable RNORMP(IPHAS).

At this level of \CS \, the array \var{TRAV} contains the second member obtained from \fort{predvv} without the user source terms. The computational procedure for the calculation of \var{RNORMP} is enumerated below:

\begin{enumerate}

\item $\displaystyle \var{TRAV}(I) = \widetilde{\vect{u}}_{\,I} - \frac {\Delta t^n_I} {\rho_I \ |\Omega_i|} \var{TRAV}(I) + \frac{(\rho_I - \rho_0 ) \Delta t^n_I }{\rho_I} \vect{g}$,
\item Call to \fort{inimas} to compute the mass flux of the vector $\var{TRAV}$ (we calculated at each face  $\displaystyle \rho_{\,ij} \ \var{TRAV}_{\,ij}\,.\,\vect{S}_{\,ij}$, where $\vect{S}$ is the surface vector ). We set the total number of sweeps (or iterations) to $1$ (\var{NSWRP} = 1), which means that there is no reconstruction of the gradients during this run through \fort{inimas} (to save computation time). The boundary condition arrays passed into \fort{inimas} contain those of the velocity $\vect{u}^n$.
\item Call to \fort{divmas} to compute the divergence at each cell of the above mass flux, which is stored in the working are \var{W1}.
\item The mass source terms stored in the array \var{SMACEL} are added to \var{W1}.
\begin{equation}
\var{W1}(I) = \var{W1}(I) - \frac{|\Omega_i|}{\rho_I} \var{SMACEL}(I)
\label {SM1}
\end{equation}
\item Call to \fort{prodsc} ($\var{RNORMP} = \sqrt{\var{W1}.\var{W1}}$). \var{RNORMP} will be used in the stop test of the iterative pressure solver to normalize the residual (see routine \fort{gradco} for the conjugate gradient inversion method).

\end{enumerate}

\etape{Preparation for solving the system}

\begin{itemize}

\item Call to \fort{grdcel} for computation of the pressure gradient $p^n$. The result is stored in \var{TRAV}. At this level, \var{TRAV} contains $\displaystyle \frac{\partial p^n}{\partial x}$, $\displaystyle \frac{\partial p^n}{\partial y}$, $\displaystyle \frac{\partial p^n}{\partial z}$.

\item Introduction of the explicit pressure cell-based gradient $p^n$ for use with the Rhie~\&~Chow filter.

$$ \var{TRAV}(I) = \widetilde{\vect{u}}_{\,I} + \frac{\var{ARAK}}{\rho_I} \ \tens{T}^n_{\,I}\ \grad{p^n}_{\,I}$$

\var{ARAK} represents the "Arakawa" coefficient, so-named within the code, whose default setting is $1$ although this value can be modified by the user in \fort{usini1}. To simplify the notations, we define $\var{ARAK} = \alpha_{\,Arak}$.

\item Call to \fort{inimas} which calculates the mass flux in \var{TRAV}. The boundary conditions applied in this case are those of the velocity (\textit{cf.} subroutine \fort{navstv}). This is still just an approximation of the boundary conditions contained in \var{TRAV}. The mass flux (\textit{cf.} subroutine \fort{inimas} for further details concerning the calculation at the boundary faces) is thus equal to:

$$m_{\,ij} = \left[\rho \widetilde{\vect{u}}+ \displaystyle \alpha_{\,Arak}\ \tens{T}^n\ \grad(p^n)\right]_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij}$$

\item Call to \fort{itrmas} to increment the mass flux at the faces\footnote{$(\grad{p^n})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij}$ is the gradient normal to the face that is equal to  $\displaystyle \frac{ p^n_J - p^n_I}{\overline{IJ}} S_{\,ij}$ on an orthogonal mesh.} by $$- \alpha_{\,Arak}\  \widetilde{T}^{\,n}_{\,ij}\ (\grad{p}^n)_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij}.$$

\item Call to \fort{clmlga} to compute the inversion of the pressure matrix with the algebraic multigrid algorithm.

\item initialisation of $\delta p$, $\delta(\delta p)$ and \var{SMBR} to 0. \var{SMBR} will serve to store the second member. In the code, $\delta p$ and $\delta(\delta p)$ are contained respectively in \var{RTP(*,IPRIPH)} and \var{DRTP}.
\item Call to \fort{divmas} for calculation of the divergence of the mass flux resulting from the last call to \fort{itrmas}. This divergence is stored in the working array \var{W7}.
\item Addition of the contributions of the mass source terms\footnote{The array \var{W7} contains the second member without the gradient of $\delta p$. It therefore remains invariant at each sweep. On the other hand, \var{SMBR} contains the entire second member and consequently varies at each sweep.} to \var{W7}.
\begin{equation}
\var{W7}(I) =  \var{W7}(I) - |\Omega_i|\ \var{SMACEL}(I)
\label {SM2}
\end{equation}
\end{itemize}

\etape{Loop over the non orthogonalities}
Assuming that the mesh is orthogonal, a single inversion would enable resolving the problem. The loop over the non-orthogonal elements is described below:

\begin{itemize}

\item Start of the loop at $k$ (thereafter, we are at $k+1$)

\begin{itemize}

\item[$\star$] Updating of \var{SMBR} at the start of the loop\footnote{The sign "-" results from the construction of the matrix.}.
$$ \var{SMBR}(I) = -\var{W7}(I) - \var{SMBR}(I) $$

\item[$\star$] Calculation of the norm of \var{SMBR} in \fort{prodsc}. It is called \var{RESIDU} in the code. As we solve the system incrementally, the second member must cancel out convergence.
\item[$\star$] If $ \var{RESIDU} < 10 \times \varepsilon \times \var{RNORMP}$, convergence is attained\footnote{$\varepsilon$ is the pressure convergence tolerance which can be modified by the user in \fort{usini1}, {\it via} the array \var{EPSILO}.}.
\begin{itemize}
\item[$\Rightarrow$] Call to \fort{itrmas} to reupdate the mass flux with the facet gradient $(\grad (\delta p)^k)_{\,f}$. We calculate at each face $\widetilde{T}^{\,n}_{\,ij}\ (\grad(\delta p)^k)_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij}$ et $\widetilde{T}^{\,n}_{\,b_{\,ik}}\ (\grad(\delta p)^k)_{\,f_{\,b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}}$.

\item[$\Rightarrow$] Reupdate\footnote{$(\delta p)^k = \sum\limits_{l=1}^{l=k} (\delta(\delta p))^l $ is stored in \var{RTP(*,IPRIPH)}.} of the pressure $p^{n+1} = p^{n} + \sum\limits_{l=1}^{l=k} (\delta(\delta p))^l$.
\end{itemize}

\item[$\star$] Alternatively,
\begin{itemize}
\item[$\Rightarrow$] $(\delta(\delta p))^{k+1} = 0$,
\item[$\Rightarrow$] Call to \fort{invers} for inversion of the system (\ref{Base_Resopv_equation resolue}). The inversion algorithm stop test data \var{RESIDU} is normalized by \var{RNORMP} (see \fort{gradco} for the inversion of the pressure operator).
\end{itemize}
\item[$\star$] If the maximum number of sweeps is attained,
\begin{itemize}
\item[$\Rightarrow$] Call to \fort{itrmas} to increment the mass flux by the pressure gradient $(\delta p)^{k}$.
\item[$\Rightarrow$] Second call to \fort{itrmas} for incrementation of the mass flux with the non-reconstructed gradient of $(\delta(\delta p))^{k+1}$ to ensure a final divergence-free field is obtained, thereby assuring consistency with the pressure matrix which does not take non orthogonalities into account\footnote{Reference can be made to the subroutine \fort{navstv} for further detail.}.
\item[$\Rightarrow$] Update of the pressure increment $(\delta p)^{k+1} = (\delta p)^{k} + (\delta(\delta p))^{k+1}$.
\end{itemize}
\item[$\star$]  Otherwise,
\begin{itemize}
\item[$\Rightarrow$] Incrementation of the mass flux taking into account a coefficient of relaxation. $ (\delta p)^{k+1} = (\delta  p)^{k} + \var{RELAX} \times (\delta(\delta) p)^{k+1}$. The relaxation factor has a default setting of $1$; however this can be modified in \fort{usini1}.
\item[$\Rightarrow$] Call to \fort{itrgrp} for calculation of the $\tens{T}^n\ \grad{(\delta p)}$ part of the second member (to which the array \var{W7} will be added at the start of a (new) loop).
$$\var{SMBR}(I) = \sum\limits_{j\in Neigh(i)}\,\widetilde{T}^{\,n}_{\,ij}\ (\grad(\delta p)^k)_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij} + \sum\limits_{k\in {\gamma_b(i)}} \,\widetilde{T}^n_{\,b_{\,ik}}\ (\grad{(\delta p)^k})_{\,f_{b_{\,ik}}}\,.\,\vect{S}_{\,b_{\,ik}}$$
\end{itemize}
\end{itemize}

\item End of the loop
\end{itemize}
\etape {Updating of the pressure }
We update the pressure with the sum of the increments of $\delta p$.
$$p^{n+1} = p^n + (\delta p)_{k_{conv}}$$
where,
$$(\delta p)_{k_{conv}} = \sum\limits_{k=1}^{k=k_{conv}}{(\delta(\delta p))^k} $$

In practice, \var{RTP} contains $(\delta p)^k_{conv}$. We therefore increment $p^n$ by \var{RTP}:

$$\var{RTP} = \var{RTPA} + \var{RTP} $$

\clearpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points to treat}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
There are a number of outstanding issues that still need to be resolved:
\begin{enumerate}

\item The use of the normal gradient as an approximation of the pressure increment gradient can pose problems in terms of consistency, as indicated notably in the remark below.

\item The approximation $\ \widetilde{T}^n \approx (\tens{T}^n \ \vect{n})\,.\,\vect{n}\ $ is not made for the reconstruction of the gradients in the second member of the pressure equation. Nor is it made when calculating the cell-based gradient during application of the Rhie \& Chow filter.

\item Use of the weighting factor $\displaystyle \frac {1}{2}$ to improve numerical stability when computing calculations based on the values at the faces  (e.g. in \fort{itrmas} or \fort{itrgrp} during reconstruction of the pressure increment gradient at the face).

\item We could verify the normalized residual calculation (see \fort{predvv}).

\item When computing the mass flux of $\displaystyle \widetilde{u}+
\frac{\alpha}{\rho}\ \tens{T}\ \grad{p^n}$ in \fort{inimas}, we use the boundary conditions of the velocity at the time level $n$. The validity of this approach remains to be clarified, particularly for the
boundary conditions at the outlet. More generally, further analysis of the boundary conditions of the variables in
\fort{navstv} is required. This issue is linked to the problem highlighted at the end of \fort{visecv}.

\item During the convergence test for the loop over the non-orthogonal elements, we multiply the tolerance by a factor of 10. Is this really necessary?

\item The problem with the relaxation factor used during updates of the pressure-correction term in the loop over the non orthogonalities (it might be worthwhile replacing this with a dynamic relaxation factor).

\item Use of the Rhie \& Chow filter can prove quite problematic in some configurations. However, before undertaking any work to modify this, it would be worthwhile to first verify its utility by assessing whether or not it plays a clear and positive role in any of the configurations.

\end{enumerate}