%                      Code_Saturne version 1.3
%                      ------------------------
%
%     This file is part of the Code_Saturne Kernel, element of the
%     Code_Saturne CFD tool.
% 
%     Copyright (C) 1998-2007 EDF S.A., France
%
%     contact: saturne-support@edf.fr
% 
%     The Code_Saturne Kernel 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.
% 
%     The Code_Saturne Kernel 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 the Code_Saturne Kernel; if not, write to the
%     Free Software Foundation, Inc.,
%     51 Franklin St, Fifth Floor,
%     Boston, MA  02110-1301  USA
%
%-----------------------------------------------------------------------
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Mise en \oe uvre}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
On expose dans ce qui suit l'algorithme tel qu'il est �crit dans \fort{resolp}.

$\vect{T}^n$ d�signera un tableau de dimension $3$ contenant les pas de temps locaux dans chaque direction (pour l'utilisation du "couplage vitesse-pression renforc\'e"). Dans le cas de l'algorithme "couplage faible vitesse-pression", on garde cette notation, les pas de temps �tant �gaux dans les trois directions de l'espace. 

\etape {Calcul de la matrice du syst�me � r�soudre}
\begin{itemize}
\item calcul du coefficient de diffusion aux faces pour le Laplacien de pression (le coefficient de diffusion utilise le pas de temps de calcul ou celui du "couplage vitesse-pression renforc\'e"). Deux cas se pr�sentent suivant l'algorithme de couplage vitesse/pression choisi par l'utilisateur :

\begin{enumerate}
\item Appel de \fort{viscfa} avec une viscosit� totale �gale au pas de temps $\Delta t_I^n$ pour l'algorithme "couplage faible vitesse-pression" (\var{IPUCOU} = 0),
\item Appel de \fort{visort} avec une viscosit� totale diagonale pour l'algorithme du "couplage vitesse-pression renforc\'e" (\var{IPUCOU} = 1). C'est \`a ce niveau qu'est calcul\'e $\widetilde{T}^n_{\,ij}$. Les pas de temps �quivalents calcul�s auparavant dans la subroutine \fort{preduv} sont stock�s dans le tableau \var{TPUCOU}.
\end{enumerate}

\item Appel de \fort{matrix} pour la construction de la matrice de diffusion de la pression (sans les termes de reconstruction qui ne peuvent pas �tre pris en compte si on veut garder une structure de matrice creuse) avec la viscosit� calcul�e pr�c\'edemment et le tableau \var{COEFB} des conditions aux limites de la pression $p^n$ (on impose une condition de Neumann homog\`ene (resp. Dirichlet homog\`ene) sur $\delta p$ pour une condition de Neumann (resp. Dirichlet) pour $p$).

\end{itemize}

\etape{Calcul du r�sidu de normalisation \var{RNORMP}}
Cette \'etape est r\'ealis\'ee directement dans \fort{preduv} \`a partir de la version
1.1.0.s, et le r\'esidu de normalisation est transmis par l'interm\'ediaire de
la variable RNORMP(IPHAS). 

Le tableau \var{TRAV} contient � ce niveau de \CS \ le second membre issu de \fort{preduv} sans les termes source utilisateur. Le calcul de \var{RNORMP} se fait de la fa�on suivante :

\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 Appel de \fort{inimas} pour calculer le flux de masse du vecteur $\var{TRAV}$ (on calcule � chaque face $\displaystyle \rho_{\,ij} \ \var{TRAV}_{\,ij}\,.\,\vect{S}_{\,ij}$, o� $\vect{S}$ est le vecteur surface). On impose le nombre de \textit{sweeps} (ou d'it\'erations) total \`a $1$ (\var{NSWRP} = 1), ce qui signifie que l'on ne recontruit pas les gradients lors de ce passage dans \fort{inimas} (pour des raisons de gain de temps de calcul). Les tableaux des conditions aux limites qui sont pass�s dans \fort{inimas} sont ceux de la vitesse $\vect{u}^n$.
\item Appel de \fort{divmas} pour calculer la divergence par cellule du flux de masse pr�c�dent. La divergence est stock�e dans le tableau de travail \var{W1}.
\item Les termes source de masse stock�s dans le tableau \var{SMACEL} sont ajout�s \`a \var{W1}.
\begin{equation}
\var{W1}(I) = \var{W1}(I) - \frac{|\Omega_i|}{\rho_I} \var{SMACEL}(I)
\label {SM1}
\end{equation}
\item Appel de \fort{prodsc} ($\var{RNORMP} = \sqrt{\var{W1}.\var{W1}}$). \var{RNORMP} servira dans le test d'arr�t du solveur de pression pour normaliser le r�sidu (voir la routine \fort{gradco} pour l'inversion par gradient conjugu�).

\end{enumerate}

\etape{pr�paration de la r�solution du syst�me}

\begin{itemize}
 
\item Appel de \fort{grdcel} pour le calcul du gradient de la pression $p^n$. Le r�sultat est stock� dans \var{TRAV}. \var{TRAV} contient donc � ce niveau $\displaystyle \frac{\partial p^n}{\partial x}$, $\displaystyle \frac{\partial p^n}{\partial y}$, $\displaystyle \frac{\partial p^n}{\partial z}$.

\item Introduction du gradient cellule de la pression explicite $p^n$ pour l'utilisation du filtre Rhie~\&~Chow.

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

\var{ARAK} repr�sente le coefficient dit "d'Arakawa" que l'utilisateur peut modifier dans \fort{usini1} et qui vaut $1$ par d�faut. Pour simplifier les notations, on pose $\var{ARAK} = \alpha_{\,Arak}$.

\item Appel de \fort{inimas} pour calculer le flux de masse de \var{TRAV}. Les conditions aux limites appliqu�es sont celles de la vitesse (voir sous-programme \fort{navsto}). Ceci reste une approximation des conditions aux limites de \var{TRAV}. Le flux de masse est donc �gal � (voir sous-programme \fort{inimas} pour plus de d�tails sur le calcul aux faces de bord) :

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

\item Appel de \fort{itrmas} pour incr�menter le flux de masse aux faces\footnote{$(\grad{p^n})_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij}$ est le gradient normal � la face �gal, sur un maillage orthogonal, � $\displaystyle \frac{ p^n_J - p^n_I}{\overline{IJ}} S_{\,ij}$.} de $$- \alpha_{\,Arak}\  \widetilde{T}^{\,n}_{\,ij}\ (\grad{p}^n)_{\,f_{\,ij}}\,.\,\vect{S}_{\,ij}.$$

\item Appel de \fort{clmlga} pour l'utilisation d'un algorithme de multigrille alg�brique pour l'inversion de la matrice de pression.

\item initialisation de $\delta p$, $\delta(\delta p)$ et \var{SMBR} � 0. \var{SMBR} servira � stocker le second membre. Dans le code, \var{RTP(*,IPRIPH)} et \var{DRTP} contiennent respectivement $\delta p$ et $\delta(\delta p)$.
\item Appel de \fort{divmas} pour le calcul de la divergence du flux de masse issu du dernier appel de \fort{itrmas}. Cette divergence est stock�e dans le tableau de travail \var{W7}.
\item Ajout des contributions des termes source de masse\footnote{Le tableau \var{W7} contient le second membre sans le gradient de $\delta p$. Il reste donc inchang\'e \`a chaque \textit{sweep}. En revanche, \var{SMBR} contient le second membre complet qui varie \`a chaque \textit{sweep}.} � \var{W7}.
\begin{equation}
\var{W7}(I) =  \var{W7}(I) - |\Omega_i|\ \var{SMACEL}(I)
\label {SM2}
\end{equation}
\end{itemize}

\etape{Boucle sur les non orthogonalit�s}
Une seule inversion permettrait de r�soudre le probl�me en supposant que le maillage soit orthogonal. On d�crit ci-dessous la boucle sur les non orthogonalit�s.

\begin{itemize}

\item d�but de la boucle sur $k$ (dans ce qui suit, on est en $k+1$)

\begin{itemize}

\item[$\star$] Actualisation de \var{SMBR} au d�but de la boucle\footnote{Le signe "-" r�sulte de la construction de la matrice.}.
$$ \var{SMBR}(I) = -\var{W7}(I) - \var{SMBR}(I) $$

\item[$\star$] Calcul de la norme de \var{SMBR} dans \fort{prodsc}. On l'appelle \var{RESIDU}. Comme on r�sout le syst�me par incr�ment, le second membre doit s'annuler � convergence.
\item[$\star$] Si $ \var{RESIDU} < 10 \times \varepsilon \times \var{RNORMP}$, la
convergence est atteinte\footnote{$\varepsilon$ est la tol�rance associ�e � la pression qui peut �tre
modifi�e par l'utilisateur dans \fort{usini1}, {\it via} le tableau \var{EPSILO}.}.
\begin{itemize}
\item[$\Rightarrow$] Appel de \fort{itrmas} pour la r�actualisation du flux de masse avec le gradient facette $(\grad (\delta p)^k)_{\,f}$. On calcule donc � chaque 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$] r�actualisation\footnote{$(\delta p)^k = \sum\limits_{l=1}^{l=k} (\delta(\delta p))^l $ est stock�e dans \var{RTP(*,IPRIPH)}.} de la pression $p^{n+1} = p^{n} + \sum\limits_{l=1}^{l=k} (\delta(\delta p))^l$.
\end{itemize}

\item[$\star$] Sinon, 
\begin{itemize}
\item[$\Rightarrow$] $(\delta(\delta p))^{k+1} = 0$,
\item[$\Rightarrow$] Appel de \fort{invers} pour l'inversion du syst�me (\ref{Base_Resolp_equation resolue}). Le test d'arr�t de l'algorithme d'inversion \var{RESIDU} est normalis� par \var{RNORMP} (voir \fort{gradco} pour l'inversion de la pression).
\end{itemize}
\item[$\star$] Si on atteint le nombre de \textit{sweeps} maximal, 
\begin{itemize}
\item[$\Rightarrow$] Appel de \fort{itrmas} pour l'incr�mentation du flux de masse par le gradient de la pression $(\delta p)^{k}$.
\item[$\Rightarrow$] Second appel de \fort{itrmas} pour l'incr�mentation du flux de masse avec le gradient de $(\delta(\delta p))^{k+1}$ non reconstruit pour assurer une divergence finale nulle, ceci pour rester coh�rent avec la matrice de pression qui ne prend pas en compte les non orthogonalit�s\footnote{On pourra se r\'ef\'erer au sous-programme \fort{navsto} pour plus de d\'etails.}. 
\item[$\Rightarrow$] mise � jour de l'incr�ment de pression $(\delta p)^{k+1} = (\delta p)^{k} + (\delta(\delta p))^{k+1}$.
\end{itemize}
\item[$\star$]  Sinon,
\begin{itemize}
\item[$\Rightarrow$] Incr�mentation du flux de masse en prenant en compte un coefficient de relaxation. $ (\delta p)^{k+1} = (\delta  p)^{k} + \var{RELAX} \times (\delta(\delta) p)^{k+1}$. Le coefficient de relaxation est mis par d�faut � $1$ et peut �tre modifi� dans \fort{usini1}.
\item[$\Rightarrow$] Appel de \fort{itrgrp} pour le calcul de la partie en $\tens{T}^n\ \grad{(\delta p)}$ du second membre (auquel sera ajout\'e le tableau \var{W7} en d\'ebut de boucle).
$$\var{SMBR}(I) = \sum\limits_{j\in Vois(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 fin de la boucle
\end{itemize}
\etape {r�actualisation de la pression}
On r�actualise la pression par la somme des incr�ments de $\delta p$.
$$p^{n+1} = p^n + (\delta p)_{k_{conv}}$$
avec,
$$(\delta p)_{k_{conv}} = \sum\limits_{k=1}^{k=k_{conv}}{(\delta(\delta p))^k} $$

En pratique, \var{RTP} contient $(\delta p)^k_{conv}$. On incr\'emente donc $p^n$ par \var{RTP} :

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

\clearpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Points \`a traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Il reste plusieurs probl\`emes \`a r\'esoudre :
\begin{enumerate}

\item L'approximation du gradient de l'incr\'ement de la pression par le gradient normal peut poser des probl\`emes de consistance, notamment avec la remarque ci-dessous.
\item L'approximation $\ \widetilde{T}^n \approx (\tens{T}^n \ \vect{n})\,.\,\vect{n}\ $ n'est pas faite lors de la reconstruction des gradients au second membre de l'\'equation de pression. Elle n'est pas faite non plus pour le calcul du gradient cellule lors de l'application du filtre Rhie \& Chow.

\item Utilisation de pond\'eration $\displaystyle \frac {1}{2}$ lors de calculs de valeurs aux faces pour des raisons de stabilit\'e num\'erique (ex. dans \fort{itrmas} ou \fort{itrgrp} lors de la reconstruction du gradient de l'incr\'ement de pression \`a la face).

\item On pourra v\'erifier le calcul du r\'esidu de normalisation (voir \fort{preduv}).

\item Lors du calcul du flux de masse de $\displaystyle \widetilde{u}+
\frac{\alpha}{\rho}\ \tens{T}\ \grad{p^n}$ dans \fort{inimas}, on utilise les
conditions aux limites de la vitesse au temps $n$. Ceci reste un point peu clair
particuli\`erement pour les conditions aux limites de sortie. Il faudra plus
g\'en\'eralement analyser les conditions aux limites des variables dans
\fort{navsto}. Ce probl\`eme est \`a relier \`a celui signal\'e pour \fort{vissec}.

\item Lors du test de convergence de la boucle sur les non orthogonalit\'es, on multiplie la tol\'erance par 10. Est-ce r�ellement n�cessaire ?

\item Il faudrait revoir le probl\`eme de relaxation lors de l'actualisation de
la  pression dans la boucle sur les non orthogonalit\'es (un coefficient de relaxation dynamique serait peut-�tre int\'eressant).

\item Le filtre de Rhie \& Chow peut s'av�rer assez g�nant dans certaines configurations. Il serait int�ressant, avant toute
travail de modification, de v�rifier
qu'il est vraiment utile, en identifiant clairement une configuration o� il joue un r�le positif certain, si une telle configuration
existe.

\end{enumerate}

