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

La variable dont il faut calculer le gradient est contenue dans le tableau
\var{PVAR}. Les conditions aux limites associ\'ees sont disponibles au travers
des tableaux \var{COEFAP} et \var{COEFBP} qui repr\'esentent respectivement les
grandeurs $A$ et $B$ utilis\'ees ci-dessus. Les trois composantes du gradient 
seront contenues, en sortie du sous-programme, dans les tableaux \var{DPDX},
\var{DPDY} et \var{DPDZ}.  


\etape{Calcul de la matrice}
Les \var{NCEL} matrices $\tens{C}_{i}$ (matrices $3\times 3$) sont 
stock\'ees dans le tableau \var{COCG}, 
(de dimension $\text{\var{NCELET}}\times 3\times 3$). Ce dernier est initialis\'e \`a z\'ero,
puis son remplissage est r\'ealis\'e dans des boucles sur les faces internes et 
les faces de bord. Les matrices $\tens{C}_{i}$ \'etant sym\'etriques, ces boucles ne
servent qu'\`a remplir la partie triangulaire sup\'erieure, le reste \'etant
compl\'et\'e \`a la fin par sym\'etrie. 

Pour \'eviter de r\'ealiser plusieurs fois les m\^emes calculs 
g\'eom\'etriques, on conserve, en sortie de sous-programme, 
dans le  tableau \var{COCG}, l'inverse des \var{NCEL} matrices $\tens{C}_{i}$.
De plus, pour les \var{NCELBR} cellules qui ont au moins une face
de bord, on conserve dans tableau \var{COCGB}, de dimension 
$\text{\var{NCELBR}}\times 3\times 3$, la contribution aux matrices $\tens{C}_{i}$ 
des termes purement g\'eom\'etriques. On pr\'ecise ces points ci-dessous. 
Notons donc d\`es \`a pr\'esent qu'il ne faut pas utiliser les tableaux 
\var{COCG} et \var{COCGB} par ailleurs comme tableaux de travail. 


\hspace*{1cm}{\bf Cellule ne poss\'edant pas de face de bord}\\
Lorsque, pour une cellule, aucune des faces n'est une face de bord du domaine,
l'expression de la matrice $\tens{C}_{i}$ ne fait intervenir que des grandeurs
g\'eom\'etriques et elle reste inchang\'ee tant que le maillage n'est
pas d\'eform\'e. Son inverse n'est donc calcul\'e qu'une seule
fois, au premier appel de \var{GRADMC} avec \var{ICCOCG=1} (l'indicateur 
\var{INICOC}, local \`a \var{GRADMC}, est positionn\'e \`a 0 
d\`es lors que ces calculs g\'eom\'etriques ont \'et\'e r\'ealis\'es une fois).
Le tableau \var{COCG} est ensuite 
r\'eutilis\'e lors des appels ult\'erieurs au sous-programme 
\var{GRADMC}.  

\hspace*{1cm}{\bf Cellule poss\'edant au moins une face de bord}\\
Lorsque l'ensemble des faces d'une cellule contient au moins une face de bord
du domaine, un terme contributeur aux matrices  $\tens{C}_{i}$ est 
sp\'ecifique \`a la variable dont on cherche 
\`a calculer le gradient, au travers du coefficient $B_{ik}$   
issu des conditions aux limites. Il s'agit de~:

\begin{equation}\label{Base_Gradmc_eq_termes_de_bord_de_C_gradmc}
\sum\limits_{k\in \gamma_b(i)}\left(\vect{n}_{b,ik}+\frac{1-B_{ik}}{||\vect{I'F}||}\,\vect{II'}\right)_{l}
                               \left(\vect{n}_{b,ik}+\frac{1-B_{ik}}{||\vect{I'F}||}\,\vect{II'}\right)_{m}
\end{equation} 
         


Au premier appel r\'ealis\'e avec \var{ICCOCG=1}, on calcule la contribution des faces
internes et on les stocke dans le tableau \var{COCGB}, qui 
sera disponible lors des appels ult\'erieurs. En effet, la
contribution des faces internes est de nature purement g\'eom\'etrique et reste
donc inchang\'ee tant que le maillage ne subit pas de d\'eformation. Elle s'\'ecrit~:
\begin{equation}\notag
\displaystyle\sum\limits_{j\in Vois(i)}\frac{1}{||\vect{IJ}||^2}(\vect{IJ})_l(\vect{IJ})_m	
\end{equation}
\`A tous les appels réalisés avec \var{ICCOCG=1}, les termes
qui d\'ependent des faces de bord (\ref{Base_Gradmc_eq_termes_de_bord_de_C_gradmc}) sont ensuite 
calcul\'es et on additionne cette contribution et \var{COCGB} qui contient celle 
des faces internes : on obtient ainsi les matrices $\tens{C}_i$ dans le tableau
\var{COCG}. Leur inverse se calcule ind\'ependamment pour chaque
cellule et on le conserve dans \var{COCG} qui sera disponible lors des appels
ult\'erieurs. 

Lorsque \var{GRADMC} a \'et\'e appel\'e une fois avec \var{ICCOCG=1}, 
des calculs peuvent \^etre \'evit\'es en
positionnant l'indicateur \var{ICCOCG} \`a 0 (si
\var{ICCOCG} est positionn\'e \`a 1, tous les calculs relatifs aux cellules
ayant au moins une face de bord sont refaits). 

\begin{itemize}
\item [-] Si \var{GRADMC} est utilis\'e pour 
calculer le gradient de la m\^eme variable (ou, plus g\'en\'eralement, d'une
variable dont les conditions aux limites conduisent aux m\^emes valeurs du 
coefficient $B_{ik}$), les matrices $\tens{C}_{i}$ sont inchang\'ees et leur
inverse est disponible dans \var{COCG} (on positionne \var{ICCOCG} \`a 0 pour
\'eviter de refaire les calculs). 
\item [-] Dans le cas contraire, les termes relatifs aux faces de 
bord (\ref{Base_Gradmc_eq_termes_de_bord_de_C_gradmc}) sont recalcul\'es et on additionne
cette contribution et \var{COCGB} qui fournit celle 
des faces internes : on obtient ainsi les matrices $\tens{C}_{i}$ dans
\var{COCG}. Il reste alors \`a inverser ces matrices. 
\end{itemize}

\hspace*{1cm}{\bf Remarque~:}\\
Pour sauvegarder les contributions g\'eom\'etriques dans \var{COCGB}, on a
recours a une boucle portant sur les \var{NCELBR} cellules  
dont au moins une face est une face de bord du domaine. Le num\'ero de ces 
cellules est donn\'e par \var{IEL = ICELBR(II)} (\var{II} variant de 1 \`a 
\var{NCELBR}). Les op\'erations r\'ealis\'ees dans cette 
boucle sont du type \var{COCGB(II,1,1) = COCG(IEL,1,1)}. La structure 
(injective) de \var{ICELBR} permet de forcer la vectorisation de la boucle.  


\etape{Inversion de la matrice}
On calcule les coefficients de la comatrice, puis l'inverse.
Pour des questions de vectorisation, la boucle sur les \var{NCEL} \'el\'ements
est remplac\'ee par une
s\'erie de boucles en vectorisation forc\'ee sur des blocs de \var{NBLOC=1024}
\'el\'ements. Le reliquat ($\var{NCEL}-E(\var{NCEL}/1024)\times 1024$) est
trait\'e apr\`es les boucles.
La matrice inverse est ensuite stock\'ee dans \var{COCG}
(toujours en utilisant sa propri\'et\'e de sym\'etrie).

\etape{Calcul du second membre et r\'esolution}
Le second membre est stock\'e dans \var{BX}, \var{BY} et \var{BZ}. Le gradient
obtenu par r\'esolution des syst\`emes locaux est stock\'e dans \var{DPDX},
\var{DPDY} et \var{DPDZ}. 


\etape{Remarque : gradient sans reconstruction} (non consistant sur maillage non
orthogonal)\\
Dans le cas o\`u l'utilisateur souhaite ne pas reconstruire le gradient ({\it
i.e.} ne pas inclure les termes de non orthogonalit\'e au calcul du gradient),
une m\'ethode sp\'ecifique est mise en \oe uvre, qui n'a pas de rapport avec la
m\'ethode de moindres carr\'es pr\'esent\'ee ci-dessus. 

Le volume de la cellule $i$ est not\'e $\Omega_i$. 
$P_{ij}$ (resp. $P_{b,ik}$) repr\'esente la valeur
estim\'ee de la variable $P$ \`a la 
face interne $ij$ (resp. \`a la face de bord $ik$) de vecteur normal associ\'e
$\vect{S}_{ij}$ (resp. $\vect{S}_{b,ik}$). Le gradient est
simplement calcul\'e en utilisant la formule suivante~: 
\begin{equation}
\begin{array}{ll}
\vect{G}_{c,i}= 
&=\displaystyle\frac{1}{\Omega_i}\left[
\sum\limits_{j\in Vois(i)}P_{ij}\vect{S}_{ij} +
\sum\limits_{k\in\gamma_b(i)}P_{b,ik}\vect{S}_{b,ik}\right]
\end{array}
\end{equation}

Les valeurs aux faces sont obtenues simplement comme suit (avec $\displaystyle\alpha_{ij}=\frac{\overline{FJ'}}{\overline{I'J'}}$)~:
\begin{equation}
\left\{\begin{array}{ll}
P_{ij}
&= \alpha_{ij}P_i+(1-\alpha_{ij})P_j\\
P_{b,ik}
&=A_{ik} +B_{ik}\,P_i
\end{array}\right.
\end{equation}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Points \`a traiter}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\etape{Vectorisation forc\'ee}
Il est peut-\^etre possible de s'affranchir du d\'ecoupage en boucles de 1024 si 
les compilateurs sont capables
d'effectuer la vectorisation sans cette aide. On note cependant que ce
d\'ecoupage en boucles de 1024 n'a pas de co\^ut CPU suppl\'ementaire, et que
le co\^ut m\'emoire associ\'e est n\'egligeable. 
Le seul inconv\'enient r\'eside dans la relative complexit\'e  de l'\'ecriture.

\etape{Choix du vecteur $d$}
Le choix $\vect{d}_{ij}= \frac{\vect{IJ}}{||\vect{IJ}||}$ permet de calculer
simplement une composante du gradient \`a la face en ne faisant intervenir que
les valeurs de la variable au centre des cellules voisines. Le choix 
$\vect{d}_{ij}= \frac{\vect{I'J'}}{||\vect{I'J'}||}$ serait \'egalement
possible, et peut-\^etre meilleur, mais conduirait naturellement \`a faire
intervenir, pour le calcul de la composante du gradient normale aux faces, les
valeurs de la variable aux points $I'$ et $J'$, et donc les valeurs du gradient
dans les cellules voisines. Il en r\'esulterait donc un syst\`eme coupl\'e,
auquel un algorithme it\'eratif (voir \var{GRADRC}) pourrait \^etre appliqu\'e. 
L'aspect temps calcul, atout majeur de la m\'ethode actuelle, 
s'en ressentirait sans doute. 

\etape{Am\'elioration de la m\'ethode}
Cette m\'ethode rencontre des difficult\'es dans le cas de maillages assez ``non
orthogonaux'' (cas de la voiture maill\'e en t\'etra\`edres par exemple). Une
voie d'am\'elioration possible est d'utiliser un support \'etendu (le support est 
l'ensemble des cellules utilis\'ees pour calculer le gradient en une cellule
donn\'ee). Un exemple est fourni sur la figure \ref{Base_Gradmc_fig_support_gradmc} ci-dessous~: si la cellule $I$
est la cellule courante, on choisit pour
support les cellules de centre $J$ telles que la droite $(IJ)$ soit la plus
orthogonale possible \`a une face de la cellule $I$. 


\begin{figure}[htp]
\centerline{\includegraphics[height=8cm]{../Base/Gradmc/Images/support.pdf}}
\caption{\label{Base_Gradmc_fig_support_gradmc}Diff\'erents supports pour le calcul du gradient.}
\end{figure}



