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

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

\vspace{1cm}
%-------------------------------------------------------------------------------
\section*{Fonction}
%-------------------------------------------------------------------------------
Ce sous-programme, appel\'{e} entre autre par \fort{predvv}, \fort{turbke}, \fort{covofi},
\fort{resrij}, \fort{reseps}, ..., r\'{e}sout les \'{e}quations de convection-diffusion
d'un scalaire $a$ avec termes sources du 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{convection implicite}}
-\theta \ \underbrace{\dive(\mu_{\,tot}\,\grad a^{n+1})}_{\text{diffusion implicite}}
\\\\
= f_s^{\,exp}-(1-\theta) \ \underbrace{\dive((\rho \underline{u})\,a^{n})}_{\text{convection explicite}}
 + (1-\theta) \ \underbrace{\dive(\mu_{\,tot}\,\grad a^{n})}_{\text{diffusion explicite}}
\end{array}
\end{equation}
o\`{u} $\rho \underline{u}$, $f_s^{exp}$ et $f_s^{imp}$ d\'{e}signent respectivement le flux de masse, les termes sources explicites et les termes lin\'{e}aris\'{e}s en $a^{n+1}$.
$a$ est un scalaire d\'{e}fini sur toutes les cellules\footnote{$a$, sous forme discr\`ete en espace, correspond \`a un vecteur dimensionn\'e \`a \var{NCELET} de composante $a_I$, I d\'ecrivant l'ensemble des cellules.}.
Par souci de clart\'{e} on suppose, en l'absence d'indication, les propri\'{e}tes
physiques $\Phi$ (viscosit\'{e} totale $\mu_{tot}$,...) et le flux de masse $(\rho
\underline{u})$ pris respectivement aux instants $n+\theta_\Phi$ et
$n+\theta_F$, o\`{u} $\theta_\Phi$ et $\theta_F$ d\'{e}pendent des sch\'{e}mas en temps
sp\'{e}cifiquement utilis\'{e}s pour ces grandeurs\footnote{cf. \fort{introd}}.
\\
L'\'{e}criture des termes de convection et diffusion en maillage non orthogonal
engendre des difficult\'{e}s (termes de reconstruction et test de pente) qui sont
contourn\'{e}es en utilisant une m\'ethode it\'erative dont la limite, si elle
existe, est la solution de l'\'{e}quation pr\'{e}c\'{e}dente.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Discr\'etisation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Afin d'expliquer la proc\'{e}dure utilis\'{e}e pour traiter les difficult\'{e}s dues aux
termes de reconstruction et de test de pente dans les termes de
convection-diffusion, on note, de fa\c con analogue \`a ce qui est d\'efini dans
\fort{navstv} mais sans discr\'etisation spatiale associ\'ee, $\mathcal{E}_{n}$ l'op\'erateur :
\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}
L'\'equation (\ref{Base_Codits_eq_ref}) s'\'ecrit donc :
\begin{equation}
\mathcal{E}_{n}(a^{n+1}) = 0
\end{equation}
La quantit\'e  $\mathcal{E}_{n}(a^{n+1})$ comprend donc :\\
\hspace*{1.cm} $\rightsquigarrow$ $f_s^{\,imp}\,a^{n+1}$, contribution des
termes diff\'erentiels d'ordre $0$ lin\'eaire en $a^{n+1}$,\\
\hspace*{1.cm} $\rightsquigarrow$ $\theta\,\,
\dive((\rho\underline{u})\,a^{n+1})
- \theta\,\, \dive(\mu_{\,tot}\,\grad a^{n+1})$, termes de convection-diffusion
implicites complets (termes non reconstruits + termes de reconstruction)
lin\'eaires\footnote{Lors de la discr\'{e}tisation en espace, le caract\`{e}re lin\'{e}aire
de ces termes pourra cependant �tre perdu, notamment \`{a} cause du test de pente.}
en $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)$ l'ensemble
des termes explicites (y compris la partie explicite provenant du sch\'ema en
temps appliqu\'e \`a la convection diffusion).\\\\

De m\^eme, on introduit un op\'erateur $\mathcal{EM}_{n}$ approch\'e de
$\mathcal{E}_{n}$, lin\'eaire et simplement inversible, tel que son
expression contient :\\
\hspace*{1.cm}$\rightsquigarrow$ la prise en compte des termes lin\'eaires en $a$,\\
\hspace*{1.cm}$\rightsquigarrow$ la convection int\'{e}gr\'{e}e par un sch\'{e}ma d\'{e}centr\'{e} amont
(upwind) du premier ordre en espace,\\
\hspace*{1.cm}$\rightsquigarrow$ les flux diffusifs non reconstruits.\\
\begin{equation}
\mathcal{EM}_{n}(a) = f_s^{\,imp}\,a + \theta\,\,[\dive((\rho
\underline{u})\,a)]^{\textit{amont}} - \theta\,\, [\dive(\mu_{\,tot}\,\grad a)]^{\textit{N Rec}}
\end{equation}
Cet op\'erateur permet donc de contourner la difficult\'e induite par la pr\'esence d'\'eventuelles non lin\'earit\'es introduites par l'activation du test de pente lors du sch\'ema convectif, et par le remplissage important de la structure de la matrice d\'ecoulant de la pr\'esence des gradients propres \`a la reconstruction.\\
On a la relation\footnote{On pourra se reporter au sous-programme
\fort{matrix} pour plus de d\'etails relativement \`a
$\mathcal{EM_{\it{disc}}}$, op\'erateur discret agissant sur un scalaire $a$.}, pour toute cellule $\Omega_I$ de centre $I$  :
\begin{equation}\notag
\mathcal{EM_{\it{disc}}}(a,I) = \int_{\Omega_i}\mathcal{EM}_{n}(a)  \, d\Omega
\end{equation}
On cherche \`{a} r\'{e}soudre :
\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}
Soit :
\begin{equation}
\mathcal{EM}_{n}(a^{n+1}) =  \mathcal{EM}_{n}(a^{n+1}) -  \mathcal{E}_{n}(a^{n+1})
\end{equation}
On va pour cela utiliser un algorithme de type point fixe en d\'{e}finissant la
suite $(a^{n+1,\,k})_{k\in \mathbb{N}}$\footnote{Dans le cas ou le point fixe en
vitesse-pression est utilis\'{e} (\var{NTERUP}$>$ 1) $a^{n+1,0}$ est initialis\'{e} par
la derni\`{e}re valeur obtenue de $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}
o\`{u} $\delta a^{n+1,\,k+1}$ est solution de :
\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}
Soit encore, par lin\'{e}arit\'{e} de $\mathcal{EM}_{n}$ :
\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}

Cette suite, coupl\'ee avec le choix de l'op\'erateur $\mathcal{E}_{n}$, permet donc de lever la difficult\'{e} induite par la
pr\'esence de la convection (discr\'etis\'ee \`a l'aide de sch\'emas num\'eriques
qui peuvent introduire des non lin\'earit\'es) et les termes de
reconstruction. Le sch\'ema r\'eellement choisi par l'utilisateur pour la
convection (donc \'eventuellement non lin\'eaire si le test de pente est activ\'e) ainsi que les termes de
reconstruction vont \^etre pris \`a l'it\'{e}ration $k$ et trait\'es au second membre {\it via} le sous-programme \fort{bilsc2},  alors que les termes
non reconstruits sont pris \`{a} l'it\'{e}ration $k+1$ et repr\'esentent donc les
inconnues du syst\`eme lin\'eaire r\'esolu par \fort{codits}\footnote{cf. le sous-programme
\fort{navstv}.}.\\

On suppose de plus que cette suite $(a^{n+1,\,k})_k$ converge vers la solution
$a^{n+1}$ de l'\'equation (\ref{Base_Codits_Eq_ref_small}), {\it i.e.}
$\lim\limits_{k\rightarrow\infty} \delta a^{n+1,\,k}\,=\,0$, ceci pour tout $n$ donn\'e.\\
(\ref{Base_Codits_Eq_Codits}) correspond \`a l'\'equation r\'esolue par \fort{codits}. La
matrice $\tens{EM}_{\,n}$, matrice associ\'ee \`a $\mathcal{EM}_{n}$  est
 \`a inverser, les termes non lin\'eaires sont mis au second membre mais sous forme
 explicite (indice $k$ de $a^{n+1,\,k}$) et ne posent donc plus de probl\`eme.

\minititre{Remarque 1}
La viscosit\'{e} $\mu_{\,tot}$ prise dans $\mathcal{EM}_{n}$ et dans
$\mathcal{E}_{n}$  d\'{e}pend du mod\`{e}le de turbulence utilis\'{e}. Ainsi on a
 $\mu_{\,tot}=\mu_{\,laminaire} + \mu_{\,turbulent}$
dans $\mathcal{EM}_{n}$ et dans $\mathcal{E}_{n}$ sauf lorsque l'on
utilise un mod\`{e}le $R_{ij}-\varepsilon$, auquel cas on a
$\mu_{\,tot}=\mu_{\,laminaire}$.\\
Le choix de $\mathcal{EM}_{n}$ \'{e}tant  {\it a
priori} arbitraire ($\mathcal{EM}_{n}$ doit \^etre lin\'eaire et la suite
 $(a^{n+1,\,k})_{k\in\mathbb{N}}$ doit converger pour tout $n$ donn\'e), une option des mod\`{e}les
$R_{ij}-\varepsilon$ ($\var{IRIJNU}=1$) consiste \`a for\c cer  $\mu_{\,tot}^n$
dans l'expression de $\mathcal{EM}_{n}$ \`a la
valeur $\mu_{\,laminaire}^n + \mu_{\,turbulent}^n$ lors de l'appel \`a
\fort{codits} dans le sous-programme \fort{navstv}, pour l'\'etape de
pr\'ediction de la vitesse. Ceci n'a pas de sens
physique (seul $\mu_{\,laminaire}^n$ \'{e}tant cens\'{e} intervenir), mais cela peut
dans certains cas avoir un effet stabilisateur, sans que cela modifie pour
autant les valeurs de la limite de la suite $(a^{n+1,\,k})_k$.\\

\minititre{Remarque 2}
Quand \fort{codits} est utilis\'e pour le couplage instationnaire renforc\'{e}
vitesse-pression (\var{IPUCOU}=1), on fait une seule it\'{e}ration $k$ en initialisant la suite $(a^{n+1,\,k})_{k\in\mathbb{N}}$ \`{a} z\'{e}ro. Les conditions de type Dirichlet sont
annul\'{e}es (on a $\var{INC}\,=\,0$) et le second membre est \'{e}gal \`{a} $\rho |\Omega_i|$.
Ce qui permet d'obtenir une approximation de type diagonal de
$\tens{EM}_{n}$
n\'ecessaire lors de l'\'{e}tape de correction de la vitesse\footnote{cf. le sous-programme \fort{resopv}.}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L'algorithme de ce sous-programme est le suivant :\\
- d\'{e}termination des propri\'{e}t\'{e}s de la matrice $\tens{EM}_{n}$ (sym\'etrique
si pas de convection, non sym\'etrique sinon)\\
- choix automatique de la m\'{e}thode de r\'{e}solution pour l'inverser si l'utilisateur
ne l'a pas sp\'ecifi\'e pour la variable trait\'ee. La m\'ethode de Jacobi est utilis\'ee par d\'efaut pour toute variable scalaire $a$ convect\'ee. Les m\'ethodes
disponibles sont la m\'ethode du gradient conjugu\'e, celle de Jacobi, et le
bi-gradient conjugu\'e stabilis\'e ($BICGStab$) pour les matrices non
sym\'etriques. Un pr\'econditionnement diagonal est possible et utilis\'e par
d\'efaut pour tous ces solveurs except\'{e} Jacobi.\\
- prise en compte de la p\'eriodicit\'e (translation ou rotation d'un scalaire, vecteur ou tenseur),\\
- construction de la matrice $\tens{EM}_{n}$ correspondant \`{a} l'op\'erateur
lin\'eaire $\mathcal{EM}_{n}$ par appel au sous-programme
\fort{matrix}\footnote{ On rappelle que dans \fort{matrix}, la convection  est
trait\'{e}e, quelque soit le choix de l'utilisateur, avec un sch\'{e}ma d\'{e}centr\'{e} amont d'ordre 1 en espace et qu'il n'y a pas de reconstruction
pour la diffusion. Le choix de l'utilisateur quant au
sch\'{e}ma num\'{e}rique pour la convection intervient uniquement lors de l'int\'{e}gration
des termes de convection de $\mathcal{E}_{n}$, au second membre de
(\ref{Base_Codits_Eq_Codits}) dans le sous-programme \fort{bilsc2}.}. Les termes implicites correspondant \`{a}
la partie diagonale de la matrice et donc aux contributions diff\'erentielles
d'ordre $0$ lin\'eaires en $a^{n+1}$,({\it i.e} $f_s^{imp}$), sont stock\'{e}s dans le tableau \var{ROVSDT} (r\'ealis\'e en amont du sous-programme appelant \fort{codits}).\\
- cr\'{e}ation de la hi\'{e}rarchie de maillage si on utilise le multigrille
($ \var{IMGRP}\,>0 $).\\
- appel \`{a} \fort{bilsc2} pour une \'{e}ventuelle prise en compte de la
convection-diffusion explicite lorsque $\theta \ne 0$.\\
- boucle sur le nombre d'it\'{e}rations de 1 \`a $\var{NSWRSM}$ (appel\'{e} $\var{NSWRSP}$ dans \fort{codits}).
 Les it\'{e}rations sont repr\'{e}sent\'{e}es par $k$ appel\'{e}
\var{ISWEEP} dans le code et d\'efinissent les indices de la suite $(a^{n+1,\,k})_k$
et de $(\delta a^{n+1,\,k})_k$.\\
Le second membre est scind\'e en deux parties :\\
\hspace*{1cm}{\tiny$\blacksquare$}\ un terme, affine
en  $a^{n+1,\,k-1}$, facile \`a mettre \`a jour dans le cadre de la r\'esolution
par incr\'ement, et qui s'\'ecrit :
\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$}\ les termes issus de la
convection/diffusion (avec reconstruction) calcul\'{e}e par \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}

La boucle en $k$ est alors la suivante :
\begin{itemize}
\item Calcul du second membre, hors contribution des termes de
convection-diffusion explicite $\var{SMBINI}$; le second membre complet correspondant
\`a $\mathcal{E}_{n}(a^{n+1,\,k-1})$ est, quant \`a lui, stock\'e dans le
tableau $\var{SMBRP}$, initialis\'e par $\var{SMBINI}$ et compl\'et\'e par les
termes reconstruits de convection-diffusion par appel au sous-programme
\fort{bilsc2}.\\
\`A l'it\'eration $k$, $\var{SMBINI}$ not\'e  $\var{SMBINI}^{\,k}$ vaut :\\
\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$ Avant de boucler sur $k$, un premier appel au sous-programme \fort{bilsc2} avec $\var{THETAP}=1-\theta$ permet de prendre
en compte la partie explicite des termes de convection-diffusion provenant du sch\'ema en temps.
\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}
Avant de boucler sur $k$, le second membre $\var{SMBRP}^{\,0}$ est stock\'e dans le tableau $\var{SMBINI}^{\,0}$ et sert pour l'initialisation du reste du calcul.
\begin{equation}\notag
\var{SMBINI}^{\,0} =\var{SMBRP}^{\,0}
\end{equation}
\\
$\bullet$ pour $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}
On a donc :
\begin{equation}\notag
\var{SMBINI}^{\,1}\ =\ \var{SMBINI}^{\,0} - \var{ROVSDT}\,*(\,\var{PVAR}-\,\var{PVARA})
\end{equation}
et $\var{SMBRP}^{\,1}$ est compl\'et\'e par un second appel au sous-programme \fort{bilsc2} avec $\var{THETAP}=\theta$, de mani\`ere \`a ajouter dans le second membre la partie de la convection-diffusion implicite.
\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$ pour $k = 2$,\\
de fa\c con analogue, on obtient :
\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}
Soit :
\begin{equation}\notag
\var{SMBINI}^{\,2}\ =\ \var{SMBINI}^{\,1} - \var{ROVSDT}\,*\,\var{DPVAR}^{\,1}
\end{equation}
l'appel au sous-programme \fort{bilsc2}, \'etant syst\'ematiquement fait par la suite avec $\var{THETAP}=\theta$, on obtient de m\^eme :
\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}
o\`u
\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$ pour l'it\'eration $k+1$,\\
Le tableau $\var{SMBINI}^{\,k+1}$ initialise le second membre complet
$\var{SMBRP}^{\,k+1}$ auquel vont \^etre rajout\'ees les contributions
convectives et diffusives {\it via} le sous-programme \fort{bilsc2}.\\
on a la formule :
\begin{equation}\notag
\begin{array}{ll}
\var{SMBINI}^{\,k+1}\ &= \var{SMBINI}^{\,k} - \var{ROVSDT}\,*\,\var{DPVAR}^{\,k}\\
\end{array}
\end{equation}
Puis suit le calcul et l'ajout des termes de convection-diffusion reconstruits de $-\  \mathcal{E}_{n}(a^{n+1,\,k})$, par appel au sous-programme
\fort{bilsc2}. On rappelle que la convection est prise en compte \`{a} cette \'{e}tape
par le sch\'{e}ma num\'{e}rique choisi par l'utilisateur (sch\'{e}ma d\'{e}centr\'{e} amont du
premier ordre en espace, sch\'{e}ma centr\'{e} du second ordre en espace, sch\'{e}ma
d\'{e}centr\'{e} amont du second ordre S.O.L.U. ou une
pond\'{e}ration (blending) des sch\'{e}mas dits du second ordre (centr\'{e}  ou S.O.L.U.) avec le sch\'{e}ma
amont du premier ordre, utilisation \'eventuelle d'un test de pente).\\
Cette contribution (convection-diffusion) est alors
ajout\'{e}e dans le second membre  $\var{SMBRP}^{\,k+1}$ (initialis\'e par $\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 R\'esolution du syst\`{e}me lin\'{e}aire en $\delta a^{n+1,\,k+1}$ correspondant
\`a l'\'equation (\ref{Base_Codits_Eq_Codits}) par inversion de la matrice
$\tens{EM}_{n}$, en appelant le sous programme \fort{invers}.
On calcule $a^{n+1,\,k+1}$ gr\^ace \`a la formule :
\begin{equation}\notag
a^{n+1,\,k+1} =  a^{n+1,\,k} + \delta a^{n+1,\,k+1}
\end{equation}
Soit :
\begin{equation}\notag
\var{PVAR}^{\,k+1} =  \var{PVAR}^{\,k} + \var{DPVAR}^{\,k+1}
\end{equation}

\item Traitement de la p\'eriodicit\'e et du parall\'{e}lisme.
\item Test de convergence :\\
Il porte sur la quantit\'e  $||\var{SMBRP}^{\,k+1}|| < \varepsilon
||\tens{EM}_{n}(a^{n}) + \var{SMBRP}^{\,1}|| $, o\`u $||\,.\,||$ repr\'esente la
norme euclidienne.
Si le test est v\'erifi\'e, la convergence est atteinte et on sort de la
boucle sur les it\'{e}rations. La solution recherch\'{e}e est  $a^{\,n+1} = a^{n+1,\,k+1}$.\\
Sinon, on continue d'it\'erer dans la limite des it\'{e}rations impos\'{e}es par $\var{NSWRSM}$ dans \fort{usini1}.\\
En ``continu'' ce test de convergence s'\'{e}crit aussi :
\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}
Si bien que sur maillage orthogonal avec sch\'{e}ma de convection upwind et en l'absence de terme source, la suite converge en th\'{e}orie en une unique it\'{e}ration puisque par construction~:
\begin{equation}\notag
\begin{array}{ll}
||\var{SMBRP}^{\,2}||=\,0\,& < \varepsilon ||f_s^{\,exp}||
\end{array}
\end{equation}
\end{itemize}
Fin de la boucle.
\\

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Points \`a traiter}\label{Base_Codits_section4}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\etape{Approximation $\mathcal{EM}_{n}$ de l'op\'erateur
$\mathcal{E}_{n}$}
D'autres approches visant soit \`a modifier la d\'efinition de l'approxim\'ee, prise en compte du sch\'ema centr\'e sans reconstruction par exemple,
soit \`a abandonner cette voie seraient \`a \'etudier.\\

\etape{Test de convergence}
La quantit\'e d\'efinissant le  test de convergence est \'{e}galement \`{a} revoir, \'{e}ventuellement \`{a} simplifier.

\etape{Prise en compte de $T_s^{imp}$}
Lors de la r\'{e}solution de l'\'{e}quation par \fort{codits}, le tableau \var{ROVSDT} a
deux fonctions : il sert \`{a} calculer la diagonale de la matrice (par appel de
\fort{matrix}) et il sert \`{a} mettre \`{a} jour le second membre \`{a} chaque
sous-it\'{e}ration de la r\'{e}solution en incr\'{e}ments. Or, dans le cas o\`{u} $T_s^{imp}$
est positif, on ne l'int\`{e}gre pas dans \var{ROVSDT}, afin de ne pas affaiblir la
diagonale de la matrice. De ce fait, on ne l'utilise pas pour mettre \`{a} jour le
second membre, alors que ce serait tout \`{a} fait possible. Au final, on obtient
donc un terme source trait\'{e} totalement en explicite ($T_s^{exp}+T_s^{imp}a^n$),
alors que la r\'{e}solution en incr\'{e}ments nous permettrait justement de l'impliciter
quasiment totalement ($T_s^{exp}+T_s^{imp}a^{n+1,k_{fin}-1}$, o\`{u} $k_{fin}$ est
la derni\`{e}re sous-it\'{e}ration effectu\'{e}e).\\
Pour ce faire, il faudrait d\'{e}finir deux tableaux \var{ROVSDT} dans \fort{codits}.
