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

% 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{codits}\label{ap:codits}
%

\hypertarget{codits}{}

\vspace{1cm}
%-------------------------------------------------------------------------------
\section*{Function}
%-------------------------------------------------------------------------------
This subroutine, called by \fort{predvv}, \fort{turbke}, \fort{covofi},
\fort{resrij}, \fort{reseps}, amongst others, solves the convection-diffusion equations
of a given scalar $a$ with source terms of the type:
\begin{equation}\label{Base_Codits_eq_ref}
\begin{array}{c}
\displaystyle f_s^{\,imp} (a^{n+1} - a^{n}) +
\theta \ \underbrace{\dive((\rho \underline{u})\,a^{n+1})}_{\text{implicit convection}}
-\theta \ \underbrace{\dive(\mu_{\,tot}\,\grad a^{n+1})}_{\text{implicit diffusion}}
\\\\
= f_s^{\,exp}-(1-\theta) \ \underbrace{\dive((\rho \underline{u})\,a^{n})}_{\text{explicit convection}}
 + (1-\theta) \ \underbrace{\dive(\mu_{\,tot}\,\grad a^{n})}_{\text{explicit diffusion}}
\end{array}
\end{equation}
where $\rho \underline{u}$, $f_s^{exp}$ and $f_s^{imp}$ denote respectively the mass flux, the explicit source terms and the terms linearized in $a^{n+1}$.
$a$ is a scalar defined on all cells\footnote{$a$, in spatially discrete form, corresponds to a vector sized to \var{NCELET} of component $a_I$, I describing all of the cells.}.
For clarification, unless stated otherwise it is assumed that the physical properties $\Phi$ (total
viscosity $\mu_{tot}$,...) are evaluated at time $n+\theta_\Phi$ and the mass flux $(\rho\underline{u})$
at time $n+\theta_F$, with $\theta_\Phi$ and $\theta_F$ dependent on the specific time-integration schemes
selected for the respective quantities\footnote{cf. \fort{introd}}.
\\
The discretisation of the convection and diffusion terms in non-orthogonal grids creates numerical difficulties, notably with respect to the reconstruction terms and the slope test. These are circumvented by using an iterative process for which the limit, if there is one, is the solution of the previous equation.

See the \doxygenanchor{cs__equation__iterative__solve_8c.html\#codits}{programmers reference of the dedicated subroutine}
for further details.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discretisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to explain the procedure employed to address the computational
problems, owing to the use of reconstruction terms and the slope (or gradient)
test, in treating the convection-diffusion terms, we denote by $\mathcal{E}_{n}$
(similarly to the definition given in \fort{navstv} though without the associated spatial discretisation)
the operator:
\begin{equation}\label{Base_Codits_Eq_ref_small}
\begin{array}{c}
\mathcal{E}_{n}(a) = f_s^{\,imp}\,a + \theta\,\, \dive((\rho
\underline{u})\,a) - \theta\,\, \dive(\mu_{\,tot}\,\grad a)\\
- f_s^{\,exp} -  f_s^{\,imp}\,a^{n} +(1-\theta)\,\,\dive((\rho
\underline{u})\, a^n) - (1-\theta)\,\, \dive(\mu_{\,tot}\,\grad a^n)
\end{array}
\end{equation}
Hence, equation (\ref{Base_Codits_eq_ref}) is written:
\begin{equation}
\mathcal{E}_{n}(a^{n+1}) = 0
\end{equation}
The quantity  $\mathcal{E}_{n}(a^{n+1})$ thus comprises:\\
\hspace*{1.cm} $\rightsquigarrow$ $f_s^{\,imp}\,a^{n+1}$, contribution of the linear, zeroth-order
differential terms in $a^{n+1}$,\\
\hspace*{1.cm} $\rightsquigarrow$ $\theta\,\,
\dive((\rho\underline{u})\,a^{n+1})
- \theta\,\, \dive(\mu_{\,tot}\,\grad a^{n+1})$, full implicit convection-diffusion terms
(terms without flux reconstruction + reconstructed terms),
linear\footnote{During the spatial discretisation, the linearity of these terms
could however be lost, notably through the use of the slope test.}
in $a^{n+1}$,\\
\hspace*{1.cm} $\rightsquigarrow$ $f_s^{\,exp}- f_s^{\,imp}\,a^n$ et
$(1-\theta)\,\,\dive((\rho
\underline{u})\,a^n) - (1-\theta)\,\, \dive(\mu_{\,tot}\,\grad a^n)$ all of the
explicit terms (including the explicit parts obtained from the time scheme
applied to the convection diffusion equation).\\\\

Likewise, we introduce an operator $\mathcal{EM}_{n}$ that is close to
$\mathcal{E}_{n}$, linear and simple to invert, such that its
expression contains:\\
\hspace*{1.cm}$\rightsquigarrow$ the consideration of the linear terms in $a$,\\
\hspace*{1.cm}$\rightsquigarrow$ the convection integrated with a first-order upwind scheme in space,\\
\hspace*{1.cm}$\rightsquigarrow$ the diffusive flux without reconstruction.\\
\begin{equation}
\mathcal{EM}_{n}(a) = f_s^{\,imp}\,a + \theta\,\,[\dive((\rho
\underline{u})\,a)]^{\textit{upwind}} - \theta\,\, [\dive(\mu_{\,tot}\,\grad a)]^{\textit{N Rec}}
\end{equation}
This operator helps circumvent the difficulty caused by the presence of nonlinearities (which may be introduced if the slope test is activated in conjunction with the convection scheme) and the high fill-in that occurs in the matrix structure owing to the presence of reconstruction gradients.\\
We have, for each cell $\Omega_I$ of centre $I$, the following relation\footnote{For further details regarding $\mathcal{EM_{\it{disc}}}$, the discrete operator acting on a scalar $a$, the reader can refer to the subroutine
\fort{matrix}.}:
\begin{equation}\notag
\mathcal{EM_{\it{disc}}}(a,I) = \int_{\Omega_i}\mathcal{EM}_{n}(a)  \, d\Omega
\end{equation}
and want to solve the following problem:
\begin{equation}
0 =\mathcal{E}_{n}(a^{n+1}) =  \mathcal{EM}_{n}(a^{n+1}) +  \mathcal{E}_{n}(a^{n+1}) - \mathcal{EM}_{n}(a^{n+1})
\end{equation}
Hence:
\begin{equation}
\mathcal{EM}_{n}(a^{n+1}) =  \mathcal{EM}_{n}(a^{n+1}) -  \mathcal{E}_{n}(a^{n+1})
\end{equation}
In order to do so, we will use a fixed-point algorithm, defining the
sequence $(a^{n+1,\,k})_{k\in \mathbb{N}}$\footnote{If the fixed-point iterative algorithm is used
to solve the velocity-pressure system (\var{NTERUP}$>$ 1), $a^{n+1,0}$ is initialised by
the last value that was obtained for $a^{n+1}$.}:
\begin{equation}\notag
\left\{\begin{array}{l}
a^{n+1,\,0} = a^{n}\\
a^{n+1,\,k+1} = a^{n+1,\,k} + \delta a^{n+1,\,k+1}
\end{array}\right.
\end{equation}
where $\delta a^{n+1,\,k+1}$ is the solution of:
\begin{equation}
\mathcal{EM}_{n}(a^{n+1,\,k} + \delta a^{n+1,\,k+1}) = \mathcal{EM}_{n}(a^{n+1,\,k}) - \mathcal{E}_{n}(a^{n+1,\,k})
\end{equation}
Which, by linearity of $\mathcal{EM}_{n}$, simplifies out to:
\begin{equation}
\mathcal{EM}_{n}(\delta a^{n+1,\,k+1}) = - \mathcal{E}_{n}(a^{n+1,\,k})
\label{Base_Codits_Eq_Codits}
\end{equation}

By applying this sequence with the selected operator $\mathcal{E}_{n}$, we
are able to overcome the computational difficulty that arises when dealing
with the convection (discretised using numerical schemes that can lead to the
introduction of nonlinearities) and/or reconstruction terms. The user-specified scheme for the
convection (which may be nonlinear should the slope test be activated) as
well as the reconstruction terms will be evaluated at the iteration $k$ and treated on the right-hand side {\it via} the subroutine \fort{bilsc2}, while the non-reconstructed terms are taken at iteration $k+1$ and hence represent the unknowns
of the linear system solved by \fort{codits}\footnote{cf. the subroutine
\fort{navstv}.}.\\

We assume moreover that the sequence $(a^{n+1,\,k})_k$ converges to the solution
$a^{n+1}$ of equation (\ref{Base_Codits_Eq_ref_small}), {\it i.e.}
$\lim\limits_{k\rightarrow\infty} \delta a^{n+1,\,k}\,=\,0$, for any given value of $n$.\\
The equation solved by \fort{codits} is the abovementioned (\ref{Base_Codits_Eq_Codits}). The
matrix $\tens{EM}_{\,n}$, which is associated to the operator $\mathcal{EM}_{n}$, has to be
inverted; the nonlinear terms are placed on the right-hand side, but in explicit form
(index $k$ of $a^{n+1,\,k}$) and thus cease to pose a problem.

\minititre{Remark 1}
The total viscosity $\mu_{\,tot}$ considered in $\mathcal{EM}_{n}$ and in
$\mathcal{E}_{n}$ is dependent on the turbulence model used. More specifically, the viscous
term in $\mathcal{EM}_{n}$ and in $\mathcal{E}_{n}$ is, as a general rule, taken as being  $\mu_{\,tot}=\mu_{\,laminar} + \mu_{\,turbulent}$, the exception being when an $R_{ij}-\varepsilon$ model
is used, in which case $\mu_{\,tot}=\mu_{\,laminar}$.\\
The choice of $\mathcal{EM}_{n}$ being  {\it a
priori} arbitrary, ($\mathcal{EM}_{n}$ has to be linear and the sequence
 $(a^{n+1,\,k})_{k\in\mathbb{N}}$ must converge for any given $n$), one of the options in the
$R_{ij}-\varepsilon$ models ($\var{IRIJNU}=1$) consists in forcing the viscosity $\mu_{\,tot}^n$,
that appears in the expression of $\mathcal{EM}_{n}$ to take the value
$\mu_{\,laminar}^n + \mu_{\,turbulent}^n$ when \fort{codits} is called by the subroutine \fort{navstv} for the velocity prediction step.  There is no physical reason for doing so (only the laminar viscosity $\mu_{\,laminar}^n$ is supposed to appear), however this can have a stabilising effect in certain cases without there being any concomitant effect on the limit-values of the sequence $(a^{n+1,\,k})_k$.\\

\minititre{Remark 2}
When \fort{codits} is used for strengthened transient velocity-pressure coupling
(\var{IPUCOU}=1), a single iteration $k$ is made by initialising the sequence $(a^{n+1,\,k})_{k\in\mathbb{N}}$ to zero. The Dirichlet type boundary conditions are cancelled (we have $\var{INC}\,=\,0$) and the right-hand side is equal to $\rho |\Omega_i|$, which allows us to obtain a diagonal-type approximation of $\tens{EM}_{n}$, needed for the velocity correction step \footnote{cf. subroutine \fort{resopv}.}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Implementation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The algorithm for this subroutine is as follows:\\
- determination of the properties of the $\tens{EM}_{n}$ matrix (symmetric
if there is no convection, asymmetric otherwise)\\
- automatic selection of the solution method for its inversion if the user has not already
specified one for the variable being treated. The Jacobi method is used by default for every convected scalar variable $a$. The  methods available are the conjugate gradient method, Jacobi's and the
bi-conjugate gradient stabilised method ($BICGStab$) for asymmetric matrices.  Diagonal pre-conditioning can be implemented and is used by default for all of these solvers except for Jacobi.\\
- consideration of the periodicity (translation or rotation of a scalar, vector or tensor),\\
- construction of the $\tens{EM}_{n}$ matrix corresponding to the linear operator $\mathcal{EM}_{n}$ through a call to the subroutine
\fort{matrix}\footnote{ Remember that in \fort{matrix}, regardless of the user's choice, a first-order in space upwind scheme is always used to treat the convection and there is no reconstruction for
the diffusive flux. The user's choice of numerical scheme for the convection is only applied for the integration
of the convection terms of $\mathcal{E}_{n}$, on the right-hand side of
(\ref{Base_Codits_Eq_Codits}), which computed in the subroutine \fort{bilsc2}.}. The implicit terms corresponding to the diagonal part of the matrix and hence to the zeroth-order linear differential contributions in $a^{n+1}$,({\it i.e.} $f_s^{imp}$), are stored in the array \var{ROVSDT} (realized before the subroutine calls \fort{codits}).\\
- creation of the grid hierarchy if the multigrid algorithm is used
($ \var{IMGRP}\,>0 $).\\
- call to \fort{bilsc2} to take the explicit convection-diffusion into account should
 $\theta \ne 0$.\\
- loop over the number of iterations from 1 to $\var{NSWRSM}$ (which is called $\var{NSWRSP}$ in \fort{codits}).
The iterations are represented by $k$, which is designated \var{ISWEEP} in the code and defines the indices of the sequence $(a^{n+1,\,k})_k$ and of $(\delta a^{n+1,\,k})_k$.\\
The right-hand side is split into two parts:\\
\hspace*{1cm}{\tiny$\blacksquare$}\ a term that is affine in $a^{n+1,\,k-1}$, easily updated during the course of the incremental-iterative resolution procedure, which is written as:
\begin{equation}\notag
 -f_s^{\,imp} \left(\,a^{n+1,\,k-1} - a^{n+1,0}\right) + f_s^{\,exp}- (1-\theta)\,\left[\,\dive((\rho \underline{u})\,a^{n+1,0}) - \dive(\mu_{\,tot}\,\grad a^{n+1,0})\,\right]\\
\end{equation}
\\
\hspace*{1cm}{\tiny$\blacksquare$}\ the terms arising from the
convection/diffusion computation (with reconstruction) performed by \fort{bilsc2}.\\
\begin{equation}\notag
- \theta\,\left[\,\dive\left((\rho \underline{u})\,a^{n+1,\,k-1}\right)- \dive\left(\mu_{\,tot}\,\grad a^{n+1,\,k-1}\right)\right]
\end{equation}

The loop in $k$ is then the following:
\begin{itemize}
\item Computation of the right-hand side (second member) of the equation, without the contribution of
the explicit convection-diffusion terms $\var{SMBINI}$; as for the whole right-hand side corresponding
to $\mathcal{E}_{n}(a^{n+1,\,k-1})$, it is stored in the array $\var{SMBRP}$,
initialised by $\var{SMBINI}$ and completed with the reconstructed
convection-diffusion terms by a call to the subroutine \fort{bilsc2}.\\
At iteration $k$, $\var{SMBINI}$ noted  $\var{SMBINI}^{\,k}$ is equal to:\\
\begin{equation}\notag
\var{SMBINI}^{\,k}\  = f_s^{\,exp}-(1-\theta)\,\left[\,\dive((\rho \underline{u})\,a^n) - \dive(\mu_{\,tot}\,\grad a^n)\,\right]-\,f_s^{\,imp}(\,a^{n+1,\,k-1} - a^n\,) \\
\end{equation}
\\
$\bullet$ Before starting the loop over $k$, a first call to the subroutine \fort{bilsc2} with $\var{THETAP}=1-\theta$ serves to take the explicit part (from the time advancement scheme) of the convection-diffusion terms into account.
\begin{equation}\notag
\displaystyle
\var{SMBRP}^{\,0} = f_s^{\,exp} -(1-\theta)\,[\,\dive((\rho \underline{u})\,a^n) - \dive(\mu_{\,tot}\,\grad a^n)\,]\\
\end{equation}
Similarly, before looping on $k$, the right-hand side $\var{SMBRP}^{\,0}$ is stored in the array $\var{SMBINI}^{\,0}$ and serves to initialize the computation.
\begin{equation}\notag
\var{SMBINI}^{\,0} =\var{SMBRP}^{\,0}
\end{equation}
\\
$\bullet$ for $k = 1$,
\begin{equation}\notag
\begin{array}{ll}
\var{SMBINI}^{\,1}\ &=f_s^{\,exp}-(1-\theta)\,\left[\,\dive((\rho \underline{u})\,a^n) - \dive(\mu_{\,tot}\,\grad a^n)\,\right]-\,f_s^{\,imp}\,(\,a^{n+1,\,0} - a^n\,)\\
&=f_s^{\,exp}- (1-\theta)\,\left[\,\dive((\rho \underline{u})\,a^{n+1,\,0}) - \dive(\mu_{\,tot}\,\grad a^{n+1,\,0})\,\right]-f_s^{\,imp}\,\delta a^{n+1,\,0} \\
\end{array}
\end{equation}
This calculation is thus represented by a corresponding sequence of operations on the arrays:
\begin{equation}\notag
\var{SMBINI}^{\,1}\ =\ \var{SMBINI}^{\,0} - \var{ROVSDT}\,*(\,\var{PVAR}-\,\var{PVARA})
\end{equation}
and $\var{SMBRP}^{\,1}$ is completed by a second call to the subroutine \fort{bilsc2} with $\var{THETAP}=\theta$, so that the implicit part of the convection-diffusion computation is added to the right-hand side.
\begin{equation}\notag
\begin{array}{ll}
\var{SMBRP}^{\,1} & = \var{SMBINI}^{\,1} -\theta\,\left[\,\dive((\rho \underline{u})\,a^{n+1,\,0}) - \dive(\mu_{\,tot}\,\grad a^{n+1,\,0})\,\right]\\
& = f_s^{\,exp}\ - (1-\theta)\,\left[\,\dive((\rho \underline{u})\,a^{n}) - \dive(\mu_{\,tot}\,\grad a^{n})\,\right]- f_s^{\,imp}\,(a^{n+1,\,0} -a^{n}) \\
& -\theta\,\left[\,\dive((\rho \underline{u})\,a^{n+1,\,0}) - \dive(\mu_{\,tot}\,\grad a^{n+1,\,0})\,\right]\\
\end{array}
\end{equation}
$\bullet$ for $k = 2$,\\
In a similar fashion, we obtain:
\begin{equation}\notag
\begin{array}{ll}
\var{SMBINI}^{\,2}\ &=f_s^{\,exp}-(1-\theta)\,\left[\,\dive((\rho \underline{u})\,a^n) - \dive(\mu_{\,tot}\,\grad a^n)\,\right]-\,f_s^{\,imp}\,(\,a^{n+1,\,1} - a^n\,)\\
\end{array}
\end{equation}

Hence, the equivalent array formula:
\begin{equation}\notag
\var{SMBINI}^{\,2}\ =\ \var{SMBINI}^{\,1} - \var{ROVSDT}\,*\,\var{DPVAR}^{\,1}
\end{equation}
the call to the subroutine \fort{bilsc2} being systematically made thereafter with $\var{THETAP}=\theta$, we likewise obtain:
\begin{equation}\notag
\begin{array}{ll}
\var{SMBRP}^{\,2}\ &=\ \var{SMBINI}^{\,2}-\theta\left[\dive\left((\rho \underline{u})\,a^{n+1,\,1}\right)- \dive\left(\mu_{\,tot}\,\grad \,a^{n+1,\,1}\right)\right]\\
\end{array}
\end{equation}
where
\begin{equation}\notag
a^{n+1,\,1}=\var{PVAR}^{\,1}=\var{PVAR}^{\,0}+\var{DPVAR}^{\,1}=a^{n+1,\,0}+\delta a^{n+1,\,1}
\end{equation}
$\bullet$ for iteration $k+1$,\\
The array $\var{SMBINI}^{\,k+1}$ initialises the entire right side
$\var{SMBRP}^{\,k+1}$ to which will be added the convective and diffusive contributions
{\it via} the subroutine \fort{bilsc2}.\\
The array formula is given by:
\begin{equation}\notag
\begin{array}{ll}
\var{SMBINI}^{\,k+1}\ &= \var{SMBINI}^{\,k} - \var{ROVSDT}\,*\,\var{DPVAR}^{\,k}\\
\end{array}
\end{equation}
Then follows the computation and the addition of the reconstructed convection-diffusion terms of
$-\  \mathcal{E}_{n}(a^{n+1,\,k})$, by a call to the \fort{bilsc2} subroutine. Remember
that the convection is taken into account at this step by the user-specified numerical scheme
(first-order accurate in space upwind scheme, centred scheme with second-order spatial discretisation, second-order
linear upwind "SOLU" scheme or a weighted average (blending) of one of the second-order schemes (either centred or SOLU) and the first-order upwind scheme, with potential use of a slope test).\\
This contribution (convection-diffusion) is then added in to the right side of the
equation  $\var{SMBRP}^{\,k+1}$ (initialised by $\var{SMBINI}^{\,k+1}$).
\begin{equation}\notag
\begin{array}{ll}
\var{SMBRP}^{\,k+1}\ &= \var{SMBINI}^{\,k+1} - \theta\,\left[\,\dive\left((\rho \underline{u})\,a^{n+1,\,k}\right)- \dive\left(\mu_{\,tot}\,\grad a^{n+1,\,k}\right)\right]\\
& = f_s^{\,exp}-(1-\theta)\,\left[\,\dive((\rho \underline{u})\,a^n) - \dive(\mu_{\,tot}\,\grad a^n)\,\right]- f_s^{\,imp}\,(a^{n+1,\,k} -a^{n}) \\
&-\theta\,\left[\,\dive((\rho \underline{u})\,a^{n+1,k}) - \dive(\mu_{\,tot}\,\grad a^{n+1,k})\,\right]\\
\end{array}
\end{equation}

\item Resolution of the linear system in $\delta a^{n+1,\,k+1}$ corresponding
to equation (\ref{Base_Codits_Eq_Codits}) by inversion of the $\tens{EM}_{n}$ matrix, through a call to the subroutine \fort{invers}.
We calculate $a^{n+1,\,k+1}$ based on the formula:
\begin{equation}\notag
a^{n+1,\,k+1} =  a^{n+1,\,k} + \delta a^{n+1,\,k+1}
\end{equation}
Represented in array form by:
\begin{equation}\notag
\var{PVAR}^{\,k+1} =  \var{PVAR}^{\,k} + \var{DPVAR}^{\,k+1}
\end{equation}

\item Treatment of parallelism and of the periodicity.
\item Test of convergence:\\
The test involves the quantity  $||\var{SMBRP}^{\,k+1}|| < \varepsilon
||\tens{EM}_{n}(a^{n}) + \var{SMBRP}^{\,1}|| $, where $||\,.\,||$ denotes the
Euclidean norm. The solution sought is  $a^{\,n+1} = a^{n+1,\,k+1}$. If the test is satisfied, then convergence has been reached and we exit the iteration loop. \\
If not, we continue to iterate until the upper limit of iterations imposed by $\var{NSWRSM}$ in \fort{usini1} is reached.\\
The condition for convergence is also written, in continuous form, as:
\begin{equation}\notag
\begin{array}{ll}
||\var{SMBRP}^{\,k+1}||& < \varepsilon ||f_s^{\,exp}\ - \dive((\rho \underline{u})\,a^{n}) + \dive(\mu_{\,tot}\,\grad a^{n}) \\
& +[\dive((\rho \underline{u})\,a^{n})]^{\textit{amont}} + [\dive(\mu_{\,tot}\,\grad a^{n})]^{\textit{N Rec}}||\\
\end{array}
\end{equation}
As a consequence, on orthogonal mesh with an upwind convection scheme and in the absence of source terms, the sequence converges in theory in a single iteration because, by construction:
\begin{equation}\notag
\begin{array}{ll}
||\var{SMBRP}^{\,2}||=\,0\,& < \varepsilon ||f_s^{\,exp}||
\end{array}
\end{equation}
\end{itemize}
End of the loop.
\\

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points to treat}\label{Base_Codits_section4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\etape{Approximation $\mathcal{EM}_{n}$ of the operator
$\mathcal{E}_{n}$}
Alternative approaches should be investigated with the aim of either modifying the current definition
of the approximate operator so that the centred scheme without reconstruction, for example, is
taken into account, or abandoning it altogether.\\

\etape{Test of convergence}
The quantity defined as criterion for the convergence test is also needs to be reviewed and potentially simplified.

\etape{Consideration of $T_s^{imp}$}
During the resolution of the equation by \fort{codits}, the \var{ROVSDT} array has
two functions: it serves to compute the diagonal elements of the matrix (by calling
\fort{matrix}) and it serves to update the right-hand side at each
sub-iteration of the incremental resolution. However, if $T_s^{imp}$ is positive,
we don't include it in \var{ROVSDT} so as not to reduce the diagonal dominance
of the matrix. In consequence, it is not used in the update of the right-hand side,
although it would certainly be feasible to include positive values in the updates.
The outcome of this is a source term that has been computed in a fully explicit
manner ($T_s^{exp}+T_s^{imp}a^n$), whereas the advantage of the incremental
approach is precisely that it allows for almost total implicitization of the solution
($T_s^{exp}+T_s^{imp}a^{n+1,k_{fin}-1}$, where $k_{fin}$ is
the final sub-iteration executed).\\
To achieve this, would require defining two \var{ROVSDT} arrays in \fort{codits}.