Le but de ce sous-programme est de calculer la vitesse au centre des cellules à partir du flux de masse aux faces, par moindres carrés. Utilisée après l’étape de correction de pression (cf. navsto ) cette méthode est une alternative à la technique de reconstruction à partir du gradient de l’incrément de pression (technique standard). Elle est activée quand l’indicateur IREVMC vaut 1 ou 2.
On rappelle que, à la fin de l’étape de correction de pression, le flux de masse aux faces vaut :
où u est la vitesse issue de l’étape de prédiction, D ij un opérateur de gradient aux faces
et RC ij le terme d’Arakawa (cf. navsto
pour une définition précise des notations).
Une première méthode, activée par IREVMC = 2, consiste à partir directement de
(ρ u) ijn+1.S ij pour calculer un+1 par moindres carrés. Son utilisation a
montré qu’elle semblait plus diffusive que la méthode standard (par exemple, dans le cas de la cavité entraînée)
et pouvait conduire à des résultats erronés sur des maillages ne comportant pas uniquement des tétraèdres
(ou des prismes à base triangulaire en “2D”) et des pavés (hexaèdres orthogonaux).
On note que, dans la méthode ci-dessus, on est parti d’une vitesse u au centre des cellules, qu’on a projetée aux faces pour obtenir le flux de masse, et qu’on ramène au centre des cellules par moindres carrés. Fort de cette constatation, une méthode alternative est disponible, activée par IREVMC = 1. Elle consiste à n’appliquer la méthode des moindres carrés qu’à la partie −D ij(Δ tn,δ Pn+θ) +RC ij du flux de masse et à rajouter directement
u (connu au centre des cellules) au résultat obtenu1. Cette méthode donne des résultats sensiblement meilleurs.
Soit une cellule Ωi, φij le flux de masse (total ou uniquement la partie en
gradient de pression) à travers la face la
séparant d’une cellule voisine Ωj et φ bik le flux de masse (total ou uniquement la partie en
gradient de pression)à travers la face de bord bik.
L’idéal serait de pouvoir trouver un vecteur vi telle que, pour toute cellule voisine Ωj on ait :
et l’équivalent aux faces de bords, i.e. :
Comme c’est généralement impossible d’obtenir les deux égalités précédentes2, on va simplement chercher à minimiser la fonction Fi :
Pour ce faire, on dérive Fi par rapport aux trois composantes du vecteur vi,
et on résout le système 3×3 local qui résulte :
avec S i matrice carrée 3×3 d’élément S ml i courant défini par :
Le flux de masse est passé par les arguments FLUMAS et FLUMAB.
• Calcul de la matrice
Les NCEL matrices 3× 3 sont stockées dans le tableau de travail
COCG,
de dimension NCELET× 3× 3. Ce dernier est d’abord mis à zéro, puis
son remplissage se fait dans des boucles sur les faces internes et les faces de
bord. La matrice étant symétrique, ces boucles ne
servent qu’à remplir la partie triangulaire supérieure, le reste étant
rempli par symétrie à la fin.
• Inversion de la matrice
On calcule les coefficients de la comatrice, puis l’inverse.
Pour des questions de vectorisation, la boucle sur les NCEL éléments
est remplacée par une
série de boucles en vectorisation forcée sur des blocs de NBLOC=1024
éléments. Le reliquat (NCEL−E(NCEL/1024)× 1024) est
traité après les boucles.
À la fin, la matrice inverse est stockée dans COCG
(toujours en utilisant sa propriété de symétrie).
• Calcul du second membre et résolution
Le second membre est stocké dans BX, BY et BZ. La vitesse
finale est stockée dans UX, UY et UZ.
• Vectorisation forcée
Le découpage en boucles de 1024 est-il vraiment nécessaire ? Les machines
vectorielles et les compilateurs sont-ils aujourd’hui capables
d’effectuer la vectorisation sans cette aide ? On note cependant que ce
découpage en boucles de 1024 n’a pas de coût CPU supplémentaire, et un
coût mémoire négligeable. Le seul inconvénient réside dans la
complexité de l’écriture.
• Suppression de la méthode IREVMC = 2
Sur un maillage “1D” d’hexaèdres tous orthogonaux sauf une face, on peut montrer que la méthode fait apparaître
une composante de vitesse aberrante non nulle et directement déterminée par l’angle de non orthogonalité de la
face (non consistance). On pourrait donc songer à supprimer purement cette méthode, dans la mesure où elle n’est
a priori consistante que sur une classe réduite de maillages.