1%------------------------------------------------------------------------------- 2 3% This file is part of Code_Saturne, a general-purpose CFD tool. 4% 5% Copyright (C) 1998-2021 EDF S.A. 6% 7% This program is free software; you can redistribute it and/or modify it under 8% the terms of the GNU General Public License as published by the Free Software 9% Foundation; either version 2 of the License, or (at your option) any later 10% version. 11% 12% This program is distributed in the hope that it will be useful, but WITHOUT 13% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14% FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15% details. 16% 17% You should have received a copy of the GNU General Public License along with 18% this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19% Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21%------------------------------------------------------------------------------- 22 23\programme{codits}\label{ap:codits} 24% 25 26\hypertarget{codits}{} 27 28\vspace{1cm} 29%------------------------------------------------------------------------------- 30\section*{Function} 31%------------------------------------------------------------------------------- 32This subroutine, called by \fort{predvv}, \fort{turbke}, \fort{covofi}, 33\fort{resrij}, \fort{reseps}, amongst others, solves the convection-diffusion equations 34of a given scalar $a$ with source terms of the type: 35\begin{equation}\label{Base_Codits_eq_ref} 36\begin{array}{c} 37\displaystyle f_s^{\,imp} (a^{n+1} - a^{n}) + 38\theta \ \underbrace{\dive((\rho \underline{u})\,a^{n+1})}_{\text{implicit convection}} 39-\theta \ \underbrace{\dive(\mu_{\,tot}\,\grad a^{n+1})}_{\text{implicit diffusion}} 40\\\\ 41= f_s^{\,exp}-(1-\theta) \ \underbrace{\dive((\rho \underline{u})\,a^{n})}_{\text{explicit convection}} 42 + (1-\theta) \ \underbrace{\dive(\mu_{\,tot}\,\grad a^{n})}_{\text{explicit diffusion}} 43\end{array} 44\end{equation} 45where $\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}$. 46$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.}. 47For clarification, unless stated otherwise it is assumed that the physical properties $\Phi$ (total 48viscosity $\mu_{tot}$,...) are evaluated at time $n+\theta_\Phi$ and the mass flux $(\rho\underline{u})$ 49at time $n+\theta_F$, with $\theta_\Phi$ and $\theta_F$ dependent on the specific time-integration schemes 50selected for the respective quantities\footnote{cf. \fort{introd}}. 51\\ 52The 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. 53 54See the \doxygenanchor{cs__equation__iterative__solve_8c.html\#codits}{programmers reference of the dedicated subroutine} 55for further details. 56 57%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 58%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 59\section*{Discretisation} 60%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 61%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 62In order to explain the procedure employed to address the computational 63problems, owing to the use of reconstruction terms and the slope (or gradient) 64test, in treating the convection-diffusion terms, we denote by $\mathcal{E}_{n}$ 65(similarly to the definition given in \fort{navstv} though without the associated spatial discretisation) 66the operator: 67\begin{equation}\label{Base_Codits_Eq_ref_small} 68\begin{array}{c} 69\mathcal{E}_{n}(a) = f_s^{\,imp}\,a + \theta\,\, \dive((\rho 70\underline{u})\,a) - \theta\,\, \dive(\mu_{\,tot}\,\grad a)\\ 71- f_s^{\,exp} - f_s^{\,imp}\,a^{n} +(1-\theta)\,\,\dive((\rho 72\underline{u})\, a^n) - (1-\theta)\,\, \dive(\mu_{\,tot}\,\grad a^n) 73\end{array} 74\end{equation} 75Hence, equation (\ref{Base_Codits_eq_ref}) is written: 76\begin{equation} 77\mathcal{E}_{n}(a^{n+1}) = 0 78\end{equation} 79The quantity $\mathcal{E}_{n}(a^{n+1})$ thus comprises:\\ 80\hspace*{1.cm} $\rightsquigarrow$ $f_s^{\,imp}\,a^{n+1}$, contribution of the linear, zeroth-order 81differential terms in $a^{n+1}$,\\ 82\hspace*{1.cm} $\rightsquigarrow$ $\theta\,\, 83\dive((\rho\underline{u})\,a^{n+1}) 84- \theta\,\, \dive(\mu_{\,tot}\,\grad a^{n+1})$, full implicit convection-diffusion terms 85(terms without flux reconstruction + reconstructed terms), 86linear\footnote{During the spatial discretisation, the linearity of these terms 87could however be lost, notably through the use of the slope test.} 88in $a^{n+1}$,\\ 89\hspace*{1.cm} $\rightsquigarrow$ $f_s^{\,exp}- f_s^{\,imp}\,a^n$ et 90$(1-\theta)\,\,\dive((\rho 91\underline{u})\,a^n) - (1-\theta)\,\, \dive(\mu_{\,tot}\,\grad a^n)$ all of the 92explicit terms (including the explicit parts obtained from the time scheme 93applied to the convection diffusion equation).\\\\ 94 95Likewise, we introduce an operator $\mathcal{EM}_{n}$ that is close to 96$\mathcal{E}_{n}$, linear and simple to invert, such that its 97expression contains:\\ 98\hspace*{1.cm}$\rightsquigarrow$ the consideration of the linear terms in $a$,\\ 99\hspace*{1.cm}$\rightsquigarrow$ the convection integrated with a first-order upwind scheme in space,\\ 100\hspace*{1.cm}$\rightsquigarrow$ the diffusive flux without reconstruction.\\ 101\begin{equation} 102\mathcal{EM}_{n}(a) = f_s^{\,imp}\,a + \theta\,\,[\dive((\rho 103\underline{u})\,a)]^{\textit{upwind}} - \theta\,\, [\dive(\mu_{\,tot}\,\grad a)]^{\textit{N Rec}} 104\end{equation} 105This 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.\\ 106We 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 107\fort{matrix}.}: 108\begin{equation}\notag 109\mathcal{EM_{\it{disc}}}(a,I) = \int_{\Omega_i}\mathcal{EM}_{n}(a) \, d\Omega 110\end{equation} 111and want to solve the following problem: 112\begin{equation} 1130 =\mathcal{E}_{n}(a^{n+1}) = \mathcal{EM}_{n}(a^{n+1}) + \mathcal{E}_{n}(a^{n+1}) - \mathcal{EM}_{n}(a^{n+1}) 114\end{equation} 115Hence: 116\begin{equation} 117\mathcal{EM}_{n}(a^{n+1}) = \mathcal{EM}_{n}(a^{n+1}) - \mathcal{E}_{n}(a^{n+1}) 118\end{equation} 119In order to do so, we will use a fixed-point algorithm, defining the 120sequence $(a^{n+1,\,k})_{k\in \mathbb{N}}$\footnote{If the fixed-point iterative algorithm is used 121to solve the velocity-pressure system (\var{NTERUP}$>$ 1), $a^{n+1,0}$ is initialised by 122the last value that was obtained for $a^{n+1}$.}: 123\begin{equation}\notag 124\left\{\begin{array}{l} 125a^{n+1,\,0} = a^{n}\\ 126a^{n+1,\,k+1} = a^{n+1,\,k} + \delta a^{n+1,\,k+1} 127\end{array}\right. 128\end{equation} 129where $\delta a^{n+1,\,k+1}$ is the solution of: 130\begin{equation} 131\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}) 132\end{equation} 133Which, by linearity of $\mathcal{EM}_{n}$, simplifies out to: 134\begin{equation} 135\mathcal{EM}_{n}(\delta a^{n+1,\,k+1}) = - \mathcal{E}_{n}(a^{n+1,\,k}) 136\label{Base_Codits_Eq_Codits} 137\end{equation} 138 139By applying this sequence with the selected operator $\mathcal{E}_{n}$, we 140are able to overcome the computational difficulty that arises when dealing 141with the convection (discretised using numerical schemes that can lead to the 142introduction of nonlinearities) and/or reconstruction terms. The user-specified scheme for the 143convection (which may be nonlinear should the slope test be activated) as 144well 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 145of the linear system solved by \fort{codits}\footnote{cf. the subroutine 146\fort{navstv}.}.\\ 147 148We assume moreover that the sequence $(a^{n+1,\,k})_k$ converges to the solution 149$a^{n+1}$ of equation (\ref{Base_Codits_Eq_ref_small}), {\it i.e.} 150$\lim\limits_{k\rightarrow\infty} \delta a^{n+1,\,k}\,=\,0$, for any given value of $n$.\\ 151The equation solved by \fort{codits} is the abovementioned (\ref{Base_Codits_Eq_Codits}). The 152matrix $\tens{EM}_{\,n}$, which is associated to the operator $\mathcal{EM}_{n}$, has to be 153inverted; the nonlinear terms are placed on the right-hand side, but in explicit form 154(index $k$ of $a^{n+1,\,k}$) and thus cease to pose a problem. 155 156\minititre{Remark 1} 157The total viscosity $\mu_{\,tot}$ considered in $\mathcal{EM}_{n}$ and in 158$\mathcal{E}_{n}$ is dependent on the turbulence model used. More specifically, the viscous 159term 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 160is used, in which case $\mu_{\,tot}=\mu_{\,laminar}$.\\ 161The choice of $\mathcal{EM}_{n}$ being {\it a 162priori} arbitrary, ($\mathcal{EM}_{n}$ has to be linear and the sequence 163 $(a^{n+1,\,k})_{k\in\mathbb{N}}$ must converge for any given $n$), one of the options in the 164$R_{ij}-\varepsilon$ models ($\var{IRIJNU}=1$) consists in forcing the viscosity $\mu_{\,tot}^n$, 165that appears in the expression of $\mathcal{EM}_{n}$ to take the value 166$\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$.\\ 167 168\minititre{Remark 2} 169When \fort{codits} is used for strengthened transient velocity-pressure coupling 170(\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}.}. 171 172%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 173%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 174\section*{Implementation} 175%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 176%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 177The algorithm for this subroutine is as follows:\\ 178- determination of the properties of the $\tens{EM}_{n}$ matrix (symmetric 179if there is no convection, asymmetric otherwise)\\ 180- automatic selection of the solution method for its inversion if the user has not already 181specified 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 182bi-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.\\ 183- consideration of the periodicity (translation or rotation of a scalar, vector or tensor),\\ 184- construction of the $\tens{EM}_{n}$ matrix corresponding to the linear operator $\mathcal{EM}_{n}$ through a call to the subroutine 185\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 186the diffusive flux. The user's choice of numerical scheme for the convection is only applied for the integration 187of the convection terms of $\mathcal{E}_{n}$, on the right-hand side of 188(\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}).\\ 189- creation of the grid hierarchy if the multigrid algorithm is used 190($ \var{IMGRP}\,>0 $).\\ 191- call to \fort{bilsc2} to take the explicit convection-diffusion into account should 192 $\theta \ne 0$.\\ 193- loop over the number of iterations from 1 to $\var{NSWRSM}$ (which is called $\var{NSWRSP}$ in \fort{codits}). 194The 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$.\\ 195The right-hand side is split into two parts:\\ 196\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: 197\begin{equation}\notag 198 -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]\\ 199\end{equation} 200\\ 201\hspace*{1cm}{\tiny$\blacksquare$}\ the terms arising from the 202convection/diffusion computation (with reconstruction) performed by \fort{bilsc2}.\\ 203\begin{equation}\notag 204- \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] 205\end{equation} 206 207The loop in $k$ is then the following: 208\begin{itemize} 209\item Computation of the right-hand side (second member) of the equation, without the contribution of 210the explicit convection-diffusion terms $\var{SMBINI}$; as for the whole right-hand side corresponding 211to $\mathcal{E}_{n}(a^{n+1,\,k-1})$, it is stored in the array $\var{SMBRP}$, 212initialised by $\var{SMBINI}$ and completed with the reconstructed 213convection-diffusion terms by a call to the subroutine \fort{bilsc2}.\\ 214At iteration $k$, $\var{SMBINI}$ noted $\var{SMBINI}^{\,k}$ is equal to:\\ 215\begin{equation}\notag 216\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\,) \\ 217\end{equation} 218\\ 219$\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. 220\begin{equation}\notag 221\displaystyle 222\var{SMBRP}^{\,0} = f_s^{\,exp} -(1-\theta)\,[\,\dive((\rho \underline{u})\,a^n) - \dive(\mu_{\,tot}\,\grad a^n)\,]\\ 223\end{equation} 224Similarly, 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. 225\begin{equation}\notag 226\var{SMBINI}^{\,0} =\var{SMBRP}^{\,0} 227\end{equation} 228\\ 229$\bullet$ for $k = 1$, 230\begin{equation}\notag 231\begin{array}{ll} 232\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\,)\\ 233&=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} \\ 234\end{array} 235\end{equation} 236This calculation is thus represented by a corresponding sequence of operations on the arrays: 237\begin{equation}\notag 238\var{SMBINI}^{\,1}\ =\ \var{SMBINI}^{\,0} - \var{ROVSDT}\,*(\,\var{PVAR}-\,\var{PVARA}) 239\end{equation} 240and $\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. 241\begin{equation}\notag 242\begin{array}{ll} 243\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]\\ 244& = 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}) \\ 245& -\theta\,\left[\,\dive((\rho \underline{u})\,a^{n+1,\,0}) - \dive(\mu_{\,tot}\,\grad a^{n+1,\,0})\,\right]\\ 246\end{array} 247\end{equation} 248$\bullet$ for $k = 2$,\\ 249In a similar fashion, we obtain: 250\begin{equation}\notag 251\begin{array}{ll} 252\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\,)\\ 253\end{array} 254\end{equation} 255 256Hence, the equivalent array formula: 257\begin{equation}\notag 258\var{SMBINI}^{\,2}\ =\ \var{SMBINI}^{\,1} - \var{ROVSDT}\,*\,\var{DPVAR}^{\,1} 259\end{equation} 260the call to the subroutine \fort{bilsc2} being systematically made thereafter with $\var{THETAP}=\theta$, we likewise obtain: 261\begin{equation}\notag 262\begin{array}{ll} 263\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]\\ 264\end{array} 265\end{equation} 266where 267\begin{equation}\notag 268a^{n+1,\,1}=\var{PVAR}^{\,1}=\var{PVAR}^{\,0}+\var{DPVAR}^{\,1}=a^{n+1,\,0}+\delta a^{n+1,\,1} 269\end{equation} 270$\bullet$ for iteration $k+1$,\\ 271The array $\var{SMBINI}^{\,k+1}$ initialises the entire right side 272$\var{SMBRP}^{\,k+1}$ to which will be added the convective and diffusive contributions 273{\it via} the subroutine \fort{bilsc2}.\\ 274The array formula is given by: 275\begin{equation}\notag 276\begin{array}{ll} 277\var{SMBINI}^{\,k+1}\ &= \var{SMBINI}^{\,k} - \var{ROVSDT}\,*\,\var{DPVAR}^{\,k}\\ 278\end{array} 279\end{equation} 280Then follows the computation and the addition of the reconstructed convection-diffusion terms of 281$-\ \mathcal{E}_{n}(a^{n+1,\,k})$, by a call to the \fort{bilsc2} subroutine. Remember 282that the convection is taken into account at this step by the user-specified numerical scheme 283(first-order accurate in space upwind scheme, centred scheme with second-order spatial discretisation, second-order 284linear 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).\\ 285This contribution (convection-diffusion) is then added in to the right side of the 286equation $\var{SMBRP}^{\,k+1}$ (initialised by $\var{SMBINI}^{\,k+1}$). 287\begin{equation}\notag 288\begin{array}{ll} 289\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]\\ 290& = 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}) \\ 291&-\theta\,\left[\,\dive((\rho \underline{u})\,a^{n+1,k}) - \dive(\mu_{\,tot}\,\grad a^{n+1,k})\,\right]\\ 292\end{array} 293\end{equation} 294 295\item Resolution of the linear system in $\delta a^{n+1,\,k+1}$ corresponding 296to equation (\ref{Base_Codits_Eq_Codits}) by inversion of the $\tens{EM}_{n}$ matrix, through a call to the subroutine \fort{invers}. 297We calculate $a^{n+1,\,k+1}$ based on the formula: 298\begin{equation}\notag 299a^{n+1,\,k+1} = a^{n+1,\,k} + \delta a^{n+1,\,k+1} 300\end{equation} 301Represented in array form by: 302\begin{equation}\notag 303\var{PVAR}^{\,k+1} = \var{PVAR}^{\,k} + \var{DPVAR}^{\,k+1} 304\end{equation} 305 306\item Treatment of parallelism and of the periodicity. 307\item Test of convergence:\\ 308The test involves the quantity $||\var{SMBRP}^{\,k+1}|| < \varepsilon 309||\tens{EM}_{n}(a^{n}) + \var{SMBRP}^{\,1}|| $, where $||\,.\,||$ denotes the 310Euclidean 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. \\ 311If not, we continue to iterate until the upper limit of iterations imposed by $\var{NSWRSM}$ in \fort{usini1} is reached.\\ 312The condition for convergence is also written, in continuous form, as: 313\begin{equation}\notag 314\begin{array}{ll} 315||\var{SMBRP}^{\,k+1}||& < \varepsilon ||f_s^{\,exp}\ - \dive((\rho \underline{u})\,a^{n}) + \dive(\mu_{\,tot}\,\grad a^{n}) \\ 316& +[\dive((\rho \underline{u})\,a^{n})]^{\textit{amont}} + [\dive(\mu_{\,tot}\,\grad a^{n})]^{\textit{N Rec}}||\\ 317\end{array} 318\end{equation} 319As 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: 320\begin{equation}\notag 321\begin{array}{ll} 322||\var{SMBRP}^{\,2}||=\,0\,& < \varepsilon ||f_s^{\,exp}|| 323\end{array} 324\end{equation} 325\end{itemize} 326End of the loop. 327\\ 328 329%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 330%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 331\section*{Points to treat}\label{Base_Codits_section4} 332%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 333%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 334\etape{Approximation $\mathcal{EM}_{n}$ of the operator 335$\mathcal{E}_{n}$} 336Alternative approaches should be investigated with the aim of either modifying the current definition 337of the approximate operator so that the centred scheme without reconstruction, for example, is 338taken into account, or abandoning it altogether.\\ 339 340\etape{Test of convergence} 341The quantity defined as criterion for the convergence test is also needs to be reviewed and potentially simplified. 342 343\etape{Consideration of $T_s^{imp}$} 344During the resolution of the equation by \fort{codits}, the \var{ROVSDT} array has 345two functions: it serves to compute the diagonal elements of the matrix (by calling 346\fort{matrix}) and it serves to update the right-hand side at each 347sub-iteration of the incremental resolution. However, if $T_s^{imp}$ is positive, 348we don't include it in \var{ROVSDT} so as not to reduce the diagonal dominance 349of the matrix. In consequence, it is not used in the update of the right-hand side, 350although it would certainly be feasible to include positive values in the updates. 351The outcome of this is a source term that has been computed in a fully explicit 352manner ($T_s^{exp}+T_s^{imp}a^n$), whereas the advantage of the incremental 353approach is precisely that it allows for almost total implicitization of the solution 354($T_s^{exp}+T_s^{imp}a^{n+1,k_{fin}-1}$, where $k_{fin}$ is 355the final sub-iteration executed).\\ 356To achieve this, would require defining two \var{ROVSDT} arrays in \fort{codits}. 357