Les préconditionneurs existant dans AeTher

Les préconditionneurs existant dans AeTher

Plusieurs préconditionneurs existent dans AeTher. Le premier, le préconditionnement diagonal par bloc d’inconnues aux nœuds, est toujours utilisé. Un autre préconditionnement, le block-SOR ou l’ILU(0), peut être utilisé en plus de celui-ci. 4.2.1 Préconditionnement bloc diagonal Le préconditionnement diagonal approche l’inverse d’une matrice par l’inverse de sa diagonale, qui est triviale à calculer. Si l’on note D la diagonale Figure 4.1 – Décomposition LU de la diagonale par bloc de A, le système préconditionné à gauche s’écrit D−1Ax = D−1b (4.1) Cela correspond à renormer les lignes de la matrice. Le préconditionnement à droite correspond à une renormalisation des colonnes. Pour normer les lignes et les colonnes de A, si les termes de D sont tous positifs, on peut préconditionner à gauche et à droite comme suit ( D− 1 2 AD− 1 2 y = D− 1 2b x = D− 1 2 y (4.2) La matrice issue de la discrétisation par la méthode des élément finis des équations de Navier-Stokes est définie par bloc, comme exposé dans la section 2.2.1. Il est alors naturel d’extraire la matrice des blocs diagonaux. L’inverse de cette matrice diagonale par blocs se calcule par une décomposition LU de ses blocs, qui est illustrée par la figure 4.1. On obtient deux matrices, notées DL et DU, telles que DLDU = D. Le préconditionnement droite-gauche diagonal par bloc s’obtient naturellement : ( DL −1ADU −1y = DL −1b x = DU −1y (4.3) Le préconditionnement diagonal par bloc est systématiquement utilisé dans le code AeTher, avant l’application d’un autre préconditionneur. Plutôt que d’appliquer le préconditionnement diagonal par bloc à chaque produit matrice-vecteur, la matrice est modifiée. Cette stratégie est possible car le préconditionnement diagonal ne change pas le squelette de la matrice, et le produit matriciel DL −1ADU −1 peut être calculé en une seule boucle sur tous les éléments de A.

Block-SOR

La méthode SOR (Successive Over Relaxation) est une méthode itérative de résolution de système linéaire par point fixe [135]. Elle peut également être mise sous la forme d’une matrice de préconditionnement. Comme indiqué dans la section précédente, le code AeTher modifie systématiquement la matrice A par un préconditionneur bloc diagonal. Sa diagonale est donc l’identité par bloc. On peut alors décomposer additivement cette matrice comme suit : A = I + L + U (4.4) Les matrices L et U sont les termes de la matrice A respectivement situés strictement au-dessous et au-dessus de sa diagonale par bloc. Dans ce cas, la matrice de préconditionnement SOR est définie [135] : MSOR = (I + ωL) (I + ωU) (4.5) où ω > 0 est appelé paramètre de relaxation. La matrice MSOR réalise une factorisation approchée de A : MSOR = (I + ωL) (I + ωU) = I + ω(L + U) + ω 2LU (4.6) Si ω ≈ 1, alors MSOR ≈ A + LU. Cette factorisation approchée est d’autant meilleure que les termes extra-diagonaux de A (qui entrent en jeu dans le produit LU) sont petits. Si cette dernière hypothèse n’est pas satisfaite, choisir ω petit permettra de diminuer l’influence du terme ω 2LU, mais alors I + ω(L + U) ne sera plus une aussi bonne approximation de A. 4.2.3 ILU(0) Le code AeTher est également doté d’un préconditionneur de type ILU(0) (pour Incomplete LU). La décomposition LU complète d’une matrice creuse est possible. Le solveur direct MUMPS utilise cette méthode [7, 9]. Néanmoins, les matrices L et U sont alors beaucoup plus remplies que la matrice originale, et leur calcul nécessite beaucoup plus d’espace mémoire et de temps qu’une décomposition incomplète, présentée ci-dessous. La décomposition ILU(0) est présentée dans l’algorithme 8. Il consiste à réaliser un pivot de Gauss pour trouver la décomposition LU de A, en ne conservant que les termes qui rentrent dans le squelette S de A. À la fin de l’algorithme, les matrices L et U sont stockées respectivement dans les parties triangulaires inférieure et supérieure de la matrice A fournie. La version bloc de l’algorithme découle immédiatement de l’algorithme 8 en considérant que aij ∈ R 5×5 (en 2D R 4×4 ) au lieu de aij ∈ R. Une factorisation LU des blocs diagonaux permet de calculer l’inversion a −1 kk aik.