# HMG, un solveur multigrille hybride géométrique-algébrique

Nicolas Tardieu, Pascal Tremblay

► To cite this version:

Nicolas Tardieu, Pascal Tremblay. HMG, un solveur multigrille hybride géométrique-algébrique. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01409588

HAL Id: hal-01409588

https://hal.science/hal-01409588v1

Submitted on 6 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# HMG, un solveur multigrille hybride géométrique-algébrique

Nicolas Tardieu¹, Pascal Tremblay²

¹ LAMSID, EDF-R&D
1 av du Gal de Gaulle, 91128 Clamart Cedex, France
nicolas.tardieu@edf.fr

² GIREF, Université Laval
Québec (Québec), Canada G1V0A6
pascal.tremblay@giref.ulaval.ca

**Résumé** — Un solveur multigrille pour la mécanique des solides est présenté. Il dispose d'une nature hybride algébrique-géométrique dans la mesure où il exploite à la fois le maillage et la matrice du système à résoudre. Il est adapté aux modélisations complexes courantes dans ce domaine incluant du mélange de modélisation et des relations linéaires entre degrés de liberté (DDL) par l'intermédiaire de multiplicateurs de Lagrange.

**Mots clés** — Solveur itératif, multigrille, multiplicateurs de Lagrange.

## 1 Introduction

Les méthodes multigrilles sont des méthodes de résolution de systèmes linéaires issues de la discrétisation d'équations aux dérivées partielles. Elles sont particulièrement adaptées aux systèmes de nature elliptique et présentent un caractère optimal au sens où leur complexité est en $O(n)$ [3] : pour résoudre un système à $n$ inconnues, elles requièrent $n$ opérations (à une constante multiplicative près).

L'objectif de cette communication est de présenter un solveur multigrille conçu pour répondre à certaines problématiques propres à la mécanique des solides. En effet, la plupart des travaux actuels se concentrent sur les performances (principalement parallèles) des méthodes multigrilles [1, 2, 6]. Peu abordent la complexité des modèles sous-jacents. Or cette complexité est très fréquente en mécanique des solides sous la forme de mélange de modélisations (poutres, coques, milieu continu) et de relations linéaires entre degrés de liberté (collage de maillages incompatibles, conditions de répétitivité cyclique, poutre 1D plongée dans un domaine 3D). C'est cette richesse de modélisation que nous considérons dans cette communication.

## 2 Proposition d'un solveur hybride géométrique-algébrique

### 2.1 Les méthodes multigrilles

Les méthodes multigrilles tirent leur force de leur nature multi-niveaux : au lieu de se limiter à une discrétisation, elles utilisent une hiérarchie de discrétisations, de la plus fine à la plus grossière, et l'utilisent pour calculer les corrections à apporter à l'approximation de la solution. Pour ce faire, elles utilisent quatre concepts :
- la hiérarchie de discrétisations où l'on note $\mathbf{M}_0$ jusqu'à $\mathbf{M}_N$ la plus grossière ;

- la restriction, qui est l'opération qui fait passer un objet porté par le maillage $\mathbf{M}_i$ sur le maillage $\mathbf{M}_{i+1}$ ; soit un vecteur $\mathbf{r}_i$, on note formellement sa restriction $\mathbf{r}_{i+1} = \mathbf{R}_i \mathbf{r}_i$ ;

- le prolongement, qui est l'opération qui fait passer un objet porté par le maillage $\mathbf{M}_{i+1}$ sur le maillage $\mathbf{M}_i$ ; soit un vecteur $\mathbf{r}_i$, on note formellement son prolongement $\mathbf{r}_i = \mathbf{P}_i \mathbf{r}_{i+1}$. Dans le cas d'un problème symétrique défini positif, on montre que $\mathbf{P}_i^T = \mathbf{R}_i$ est un choix optimal (il minimise l'erreur dans la norme en énergie) [3] ;

- le lissage, qui consiste en l'application d'une méthode de résolution peu coûteuse (Jacobi, Gauss-Seidel par exemple) permettant de limiter rapidement les composantes de hautes fréquences du résidu (composantes du résidu variant spacialement rapidement).

La matrice du niveau $i+1$ peut alors être déterminée par la formule de Galerkin [3] :

$$A_{i+1} = \mathbf{P}_i^T \cdot A_i \cdot \mathbf{P}_i \quad (1)$$

Les méthodes multigrilles se séparent en 2 familles :

- les méthodes géométriques (GMG), qui construisent explicitement les maillages des différents niveaux et qui les utilisent pour déterminer l'opérateur de restriction ;

- les méthodes algébriques (AMG), qui utilisent les entrées de la matrice du système pour choisir celles qui doivent être promues au niveau plus grossier et pour déterminer l'opérateur de restriction.

Chacune de ces familles présente des avantages et des inconvénients : plus générales mais plus difficiles à mettre en œuvre pour les GMG, plus adaptées à l'anisotropie mais plus sensibles à la nature de l'EDP pour les AMG.

Ces considérations nous amènent à proposer une approche hybride en essayant de combiner les intérêts des deux familles de méthodes :

- géométrique : on construit explicitement les maillages des différents niveaux. Cela a pour grand intérêt de faciliter la détermination de l'opérateur de restriction indépendamment de la modélisation

- algébrique : on choisit les sommets à promouvoir au prochain niveau sur la base des entrées de la matrice. On assure ainsi un bon traitement de l'anisotropie et de l'hétérogénéité.

## 2.2 Sélection des sommets

La démarche standard des multigrilles géométriques consiste à utiliser un algorithme de "Maximal Independant Set" (MIS) pour choisir les nœuds à promouvoir [1].

Dans notre approche, nous commençons par calculer la force de connexion de chaque arête [9]. Il s'agit d'un critère classique en AMG pour identifier les entrées de la matrice à conserver au niveau suivant. Nous en faisons ici une utilisation géométrique. Pour chaque arête, on extrait la sous-matrice $\mathbf{k}$ des DDL des nœuds $i$ et $j$ à ses extrémités. Cette matrice est composée de 4 sous-matrices : $\mathbf{k}_{ii}$ et $\mathbf{k}_{jj}$, termes diagonaux des DDL propres à $i$ et $j$; $\mathbf{k}_{ij}$ et $\mathbf{k}_{ji}$, termes hors-diagonaux des DDL de couplage entre $i$ et $j$. On définit alors :

$$w = \frac{\|\mathbf{k}_{ij}\| + \|\mathbf{k}_{ji}\|}{2(\|\mathbf{k}_{ii}\| + \|\mathbf{k}_{jj}\|)} \quad (2)$$

Cette information permettra de déterminer pour chaque nœud la direction dans laquelle l'erreur sera efficacement réduite par application du lisseur. C'est dans cette direction que l'on déraffinera la maillage.

Nous stockons donc pour chaque sommet, en examinant ses arêtes incidentes, la force de connexion maximale $w_{max}$ ainsi que la direction $\mathbf{e}_{max}$ de l'arête en question. Cela permet de déterminer le degré d'anisotropie du sommet $\alpha = 1 - \frac{w_{moyen}}{w_{max}}$. On peut alors utiliser l'algorithme présenté Figure 1 pour choisir les arêtes à déraffiner. On note que le tri inclut à la fois la force de connexion
◇ Pour tous les sommets

△ Pour chaque arête incidente décrite par le vecteur **e**

○ Si $\alpha > 0.3$

$$ - e_{proj} = e_{max} \cdot e $$

$$ - w = (1 - e_{proj}) * \|e\| $$

○ Sinon

$$ - w = \|e\| $$

○ Fin Si

△ Fin Pour

△ Déraffiner les 2 arêtes de plus grand w

◇ Fin Pour

<caption>Figure 1 – Sélection des arêtes</caption>

et la longueur de l'arête : il est en effet capital de déraffiner en premier les arêtes les plus courtes pour préserver la qualité des éléments [7].

## 2.3 Algorithme de déraffinement

Le caractère géométrique du solveur présenté est lié à la création explicite des maillages des niveaux grossiers. La méthode de déraffinement de maillage utilisée est celle de la contraction d'arêtes [7]. La présence d'entités topologiques de dimensions variables impose un traitement particulier.

On commence par déterminer la dimension topologique (DT) de chaque sommet comme le *min* des dimensions des éléments auxquels il appartient. On détecte ensuite les arêtes et les sommets qu'il ne faudra jamais faire disparaître : connexion entre domaines de DT différentes, extrémités d'arêtes de DT 1, "coins" de domaines de DT 3, etc...

On procède à plusieurs passes où l'on tente de contracter les arêtes du maillage qui n'ont pas été sélectionnées pour le niveau suivant sur la base de l'algorithme de MIS ou de forces de connexion. On prend soin de marquer comme "non-contractable" toutes les arêtes incidentes à un sommet avant vu une contraction d'arête lors de la passe courante (stratégie de marquage [8]). Chaque tentative de contraction doit vérifier un ensemble de règles de manière à préserver la topologie (le sommet effacé est celui de DT la plus grande, si l'arête est de DT plus petite que le domaine alors elle doit être une arête de frontière) et la géométrie du maillage initial (la contraction ne doit pas produire des éléments inversés ou de mauvaise qualité, la contraction d'une arête de frontière ne doit pas trop modifier la courbure de cette frontière).

Quand on a atteint le taux de déraffinement cible ou quand on a atteint le nombre maximal de passes de déraffinement ou quand le taux de réussite de contractions est inférieur à un certain critère, le maillage du niveau courant est produit, les opérateurs **P** et **K** sont calculés et l'on passe au niveau suivant.

## 2.4 Traitement des relations linéaires

La plupart des logiciels d'éléments finis offre des fonctionnalités impliquant des relations linéaires entre DDL dont la contrepartie est souvent l'apparition de multiplicateurs de Lagrange $\lambda$ aux côtés des inconnues physiques **u**. De fait, la matrice globale du système s'écrit alors :

$$ C = \begin{bmatrix} \mathbf{K} & \mathbf{B}^T \\ \mathbf{B} & \mathbf{0} \end{bmatrix} \qquad (3) $$

La méthode itérative choisie doit donc supporter le caractère indéfini du système ci-dessus. Mais le préconditionneur **M** tel que $M^{-1} \cdot C \approx I$ doit aussi le supporter. On propose pour ce faire d'utiliser
<figure><img src="image_2.png" /><figcaption>Figure 2 – Déraffinement anisotrope</figcaption></figure>

le préconditionneur à gauche issu de la factorisation par blocs de C [5] :

$$ M_L = \begin{bmatrix} \mathbf{K} & \mathbf{0} \\ \mathbf{B} & -\mathbf{S} \end{bmatrix} \qquad (4) $$

où $\mathbf{S} = \mathbf{B} \cdot \mathbf{K}^{-1} \cdot \mathbf{B}^T$ est le complément de Schur. L'application de $M_L$ sur le vecteur résidu $\mathbf{r} = [\mathbf{r}_u, \mathbf{r}_\lambda]^T$ est obtenu par la résolution approximative de :

$$ \begin{cases} \mathbf{K}\mathbf{u} = \mathbf{r}_u \\ \mathbf{S}\lambda = (\mathbf{B} \cdot \mathbf{u} - \mathbf{r}_\lambda) \end{cases} \qquad (5) $$

La résolution de (5-1) est réalisée de manière approximative en appliquant notre préconditionneur multigrille tandis que celle de (5-2) est réalisée en utilisant une factorisation incomplète de niveau 0 sur l'approximation $\mathbf{S} \approx \mathbf{B} \cdot \mathrm{diag}(\mathbf{K})^{-1} \cdot \mathbf{B}^T$.

## 2.5 Aspects informatiques

Informatiquement parlant, HMG se présente sous la forme d'une librairie C++. Elle dispose en outre d'une interface en C et en Fortran77, composée de 9 fonctions qui décrivent complètement le problème à résoudre. HMG a été interfacée avec le logiciel *Code_Aster* [4] avec lequel on été réalisés les test suivants.

# 3 Illustrations numériques

## 3.1 Anisotropie

Le premier exemple a pour but de montrer l'intérêt de l'approche hybride sur un exemple 2D simple. Il s'agit de résoudre sur un carré maillé uniformément le problème de thermique anisotrope suivant :

<figure><img src="image_3.png" /><figcaption>Figure 3 – Déraffinement d'une géométrie complexe : pompe RIS</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 4 – Evolution du résidu pour différents maillages de la pompe RIS</figcaption></figure>

$$ \begin{array}{l} -au_{xx} - bu_{yy} = f \quad \text{dans } \Omega \\ u = 0 \quad \text{sur } \Gamma \end{array} $$

<table><tr><td>a=b</td><td>a»b</td></tr></table>

Dans l'exemple considéré, $a = 1$ et $b = 10^{-6}$. On utilise GMRES avec 1 itération de pré et post-lissage par SOR. La figure 2 présente les maillages des niveaux initial et grossiers. Sur un maillage comportant 70 000 DDL, le nombre d'itérations pour atteindre la convergence était de 88 si on utilise le critère MIS contre 51 si on utilise le critère hybride.

## 3.2 Géométrie Complexe

Cet exemple a pour but de montrer l'efficacité de la technique de déraffinement ainsi que la performance globale du solveur. Pour ce faire, on considère le maillage d'un organe de centrale nucléaire, la pompe RIS. La Figure 3 montre le déraffinement du maillage, qui passe de 7 682 à 1 679 sommets (soit une réduction de 78%) en respectant la topologie et la géométrie sans pour autant disposer d'information de type CAD.

La Figure 4 présente l'évolution du résidu relatif préconditionné pour la résolution d'un problème d'élasticité linéaire sur plusieurs maillages comptant respectivement 23 046, 139 461 et 938 424 DDL. On utilise GMRES avec 1 itération de pré et post-lissage par SOR. Sur le maillage le plus fin, le temps CPU total de résolution était de 575s tandis qu'avec un solveur direct multi-frontal, il était de 779s.

<figure><img src="image_5.png" /><figcaption>Figure 5 – Evolution du résidu avec et sans relations linéaires</figcaption></figure>
### 3.3 Multiplicateurs de Lagrange

Cet exemple a pour but de montrer le traitement des multiplicateurs de Lagrange. On considère le cas simple d'un cube élastique encastré à sa base et soumis à une pression sur sa face supérieure. On étudie 2 cas : celui décrit précédemment et celui où tous les nœuds de la face supérieure subissent la contrainte d'avoir tous le même déplacement. Ceci change profondément le problème initial : il perd sa définie-positivité à cause de l'apparition de 720 multiplicateurs de Lagrange parmis les 4 713 DDL totaux. L'évolution des résidus relatifs préconditionnés pour ces 2 problèmes sont présentés sur la Figure 5. On constate l'effet important de l'apparition des multiplicateurs. Au bout d'une dizaine d'itérations, on observe une cassure de pente, interprétée comme une perte de vitesse de convergence liée à la pauvreté de l'approximation du complément de Schur.

## 4 Conclusions

Dans cette communication est proposé un solveur multigrille hybride géométrique-algébrique. Les exemples présentés montrent la qualité du déraffinement de maillage ainsi que la pertinence du traitement de l'anisotropie. Le traitement des relations linéaires introduites par des multiplicateurs de Lagrange mérite encore du travail de manière à améliorer l'approximation du complément de Schur.

## Références

[1] M.F. Adams - Parallel multigrid solvers for 3D unstructured finite element problems in large deformation elasticity and plasticity, *Int. J. Numer. Meth. Engng*, **48**, 2000

[2] M. Brezina, C. Tong, R. Becker - Parallel algebraic multigrids for structural mechanics, *SIAM Journal on Scientific Computing*, **27**, no5, 2006

[3] W.L. Briggs, V.E. Henson, S.F. McCormick - A multigrid tutorial, *Society for Industrial and Applied Mathematics*, 2000

[4] Code_Aster - http://www.code-aster.org

[5] A. El maliki, R. Guénette - Efficient preconditioning techniques for finite-element quadratic discretization arising from linearized incompressible Navier-Stokes, *Int. J. Numer. Meth. Fluids*, submitted

[6] A. Namazifard, I. D. Parsons - A distributed memory parallel implementation of the multi-grid method for solving three-dimensional implicit solid mechanics problems, *Int. J. Numer. Meth. Engng*, **61**, 2004

[7] X. Li, M.S. Shephard, M.W. Beall - 3D anisotropic mesh adaptation by mesh modification, *Computer Methods in Applied Mechanics and Engineering*, **194**(48-49), pp.4915-4950, 2005

[8] P. Tremblay - 2-D, 3-D and 4-D anisotropic mesh adaptation for the time-continuous space-time finite element method with applications to the incompressible Navier-Stokes equations, *Thèse de l'université d'Ottawa*, 2007

[9] P. Vanek, J. Mandel, M. Brezina - Algebraic Multigrid Based on Smoothed Aggregation for Second and Fourth Order Problems, *Computing*, **56**, pp. 179-196, 1996