# Pilotage de stratégies de calcul par décomposition de domaine par des objectifs de précision sur des quantités d'intérêt

Valentine Rey, Christian Rey, Pierre Gosselet

► To cite this version:

Valentine Rey, Christian Rey, Pierre Gosselet. Pilotage de stratégies de calcul par décomposition de domaine par des objectifs de précision sur des quantités d'intérêt. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01515052

HAL Id: hal-01515052

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

Submitted on 27 Apr 2017

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 CC0 1.0 - Universal - International License
# Pilotage de stratégies de calcul par décomposition de domaine par des objectifs de précision sur des quantités d'intérêt

V. Rey<fnref n="1" />, C. Rey<fnref n="2" />, P. Gosselet<fnref n="1" />

**Résumé** — Ce papier présente une borne stricte et garantie de l'écart entre la solution exacte d'un problème et la solution calculée dans le cadre d'une résolution élément fini par une méthode de décomposition de domaine. Cette borne sépare l'erreur de discrétisation (due au maillage) de l'erreur algébrique (due au solveur itératif). Cette séparation permet un pilotage du solveur et évite des itérations inutiles. Ce résultat est étendu à l'estimation d'erreur sur les quantités d'intérêt. Un enrichissement local est également proposé afin d'améliorer la borne obtenue.

**Mots clés** — vérification, décomposition de domaine, quantités d'intérêt, solveur de Krylov

## 1 Introduction

Les méthodes de décomposition de domaine [1] offrent un cadre efficace et robuste pour le traitement de problèmes discrétisés à plusieurs millions de degrés de liberté. Profitant de la puissance sans cesse accrue des ordinateurs, elles exploitent le parallelisme des architectures des calculateurs et possèdent les propriétés d'extensibilité numérique. La vérification est un domaine qui a connu de nombreux développements ces trente dernières années et propose des outils permettant d'estimer l'écart entre la solution exacte inconnue du problème continu $u_{ex}$ et la solution calculée. Parmi ces techniques, les estimateurs basés sur l'erreur en relation de comportement [2] permettent de fournir des bornes garanties ne nécessitant pas le calcul de constantes et sont disponibles pour une large gamme de problèmes.

Ce papier présente des techniques de pilotage de calcul par décomposition de domaine par des objectifs de précision sur des quantités d'intérêt. Il s'appuie sur un estimateur d'erreur fournissant une borne stricte et garantie de l'erreur exacte tout en séparant l'erreur algébrique (due à l'utilisation d'un solveur itératif) de l'erreur de discrétisation (due à la résolution par la méthode des éléments finis) [6]. Cet estimateur permet d'adapter le critère d'arrêt du solveur à la discrétisation, ce qui évite la sur résolution et donc les itérations inutiles. Dans [7], l'utilisation de cet estimateur est étendue aux quantités d'intérêt. Il est montré comment reconstruire les champs admissibles pour le problème de référence et le problème adjoint en se basant sur [4]. Les bornes classiques sur les quantités d'intérêt sont réécrites de manière à séparer les différentes sources d'erreur. Afin d'améliorer la borne, un enrichissement local de la discrétisation est proposé. Enfin, la discrétisation de l'interface restant inchangée via ce processus, des techniques de recyclage d'espace de Krylov sont envisagées.

Dans la suite de ce résumé, nous posons le cadre de travail (problème de référence, notations) et explicitons et illustrons les résultats et méthodes annoncés ci-dessus.

## 2 Notations et résultats

Considérons l'équilibre statique d'une structure occupant l'ouvert $\Omega \subset \mathbb{R}^d$ soumise à une force volumique $f$ dans $\Omega$, à un effort de traction $g$ sur $\partial_g \Omega$ et à un déplacement imposé $u_d$ sur la partie complémentaire du bord $\partial_u \Omega \neq \emptyset$. Nous nous plaçons dans l'hypothèse des petites perturbations et nous considérons le matériau linéaire élastique, caractérisé par le tenseur d'élasticité de Hooke $H$. $u$ est le champ de déplacement, $\varepsilon(u)$ la partie symétrique du gradient de $u$, $\sigma$ le tenseur de contraintes de Cauchy. Nous introduisons :

— Le sous-espace affine des champs cinématiquement admissibles $CA = \{u \in (H^1(\Omega))^d, u = u_d \text{ sur } \partial_u \Omega\}$

<footnote n="1">¹ LMT-Cachan, Ecole Normale Supérieure de Cachan, {valentine.rey,pierre.gosselet,}@lmt.ens-cachan.fr</footnote>

<footnote n="2">² SAFRAN, 1, rue Geneviève Aubé, CS 80112, 78772 Magny les Hameaux, France, christian.rey@lmt.ens-cachan.fr</footnote>
et notons $CA^0$ l'espace vectoriel associé.

— Le sous-espace affine des champs statiquement admissibles

$$ SA = \left\{ \underline{\tau} \in (L^2(\Omega))_{\text{sym}}^{d \times d} ; \forall \underline{v} \in CA^0, \int_{\Omega} \underline{\tau} : \underline{\underline{\varepsilon}}(\underline{v}) d\Omega = \int_{\omega} \underline{f} \cdot \underline{v} d\Omega + \int_{\partial_g \Omega} \underline{g} \cdot \underline{v} dS \right\} \quad (1) $$

— L'erreur en relation de comportement $e_{CR\Omega}(\underline{u}, \underline{\underline{\sigma}}) = \|\underline{\underline{\sigma}} - \mathbb{H} : \underline{\underline{\varepsilon}}(\underline{u})\|_{\mathbb{H}^{-1}, \Omega}$

$$ \text{avec } \|\underline{x}\|_{\mathbb{H}^{-1}, \omega} = \sqrt{\int_{\omega} (\underline{x} : \mathbb{H}^{-1} : \underline{x}) d\omega} $$

Le problème mécanique posé sur $\Omega$ s'écrit :

$$ \text{Trouver } (\underline{u}_{ex}, \underline{\underline{\sigma}}_{ex}) \in CA \times SA \text{ tel que } e_{CR\Omega}(\underline{u}_{ex}, \underline{\underline{\sigma}}_{ex}) = 0 $$

La solution de ce problème, dite solution exacte, existe et est unique. Au domaine $\Omega$, nous associons le sous-espace de dimension finie $CA_H$ de CA et résolvons le problème élément fini.

La mesure de la qualité de la solution via l'erreur en relation de comportement consiste à construire des champs admissibles à partir des champs obtenus par la méthode des éléments finis et mesurer la non-vérification de la relation de comportement. Puisque $CA_H \subset CA$, le champ $\underline{u}_H$, solution du problème éléments finis, est cinématiquement admissible. La construction de champ de contrainte statiquement admissible (respectant l'équilibre global ainsi que l'équilibre des efforts sur chaque arête) est plus complexe. Néanmoins, différentes techniques sont disponibles [5]. Une fois ce champ statiquement admissible $\hat{\underline{\sigma}}$ reconstruit, on obtient une borne supérieure de l'erreur $\|\underline{u}_{ex} - \underline{u}_H\|_{\mathbb{H}^{-1}, \Omega} \le e_{CR\Omega}(\underline{u}_H, \hat{\underline{\sigma}})$.

Les méthodes de décomposition de domaine sans recouvrement consistent à découper la structure en sous-domaines $\Omega^{(s)}$ formant une partition de $\Omega$. Le problème élément fini est reformulé en terme de quantités d'interface :

— Dans l'apporche primale, une unique inconnue d'interface de type déplacement est introduite si bien que la continuité des déplacements est toujours vérifiée. A chaque itération du solveur, le saut d'effort aux interfaces est réduit jusqu'à convergence (algorithme BDD [10])

— Dans l'approche duale, une unique inconnue d'interface de type effort (réaction nodale) est introduite si bien que l'équilibre de l'interface est toujours vérifié. A chaque itération du solveur, le saut de déplacement aux interfaces est réduit jusqu'à convergence (algorithme FETI [9])

Dans [4], pour ces deux algorithmes, il a été montré comment reconstruire, à chaque itération, les champs suivants :

— $\underline{u}_D = (\underline{u}_D^{(s)})_s \in CA(\Omega)$ : champ de déplacement continu sur toute la structure

— $\underline{u}_N = (\underline{u}_N^{(s)})_s$ : champ de déplacement continu par sous-domaine mais pas aux interfaces

— $\lambda_N^{(s)}$ : réactions nodales associées à $\underline{u}_N^{(s)}$ et équilibrées aux interfaces.

— $\underline{\underline{\sigma}}_N^{(s)} = \mathbb{H} : \underline{\underline{\varepsilon}}(\underline{u}_N^{(s)})$ : champ de contrainte qui, combiné à $\lambda_N^{(s)}$, permet de construire un champ statiquement admissible $\hat{\underline{\sigma}}_N^{(s)}$ en utilisant les méthodes dédiées.

Ainsi, quel que soit le niveau de convergence du solveur, il est possible de construire un couple $(\underline{u}_D, \hat{\underline{\sigma}}_N)$ de champs respectivement cinématiquement et statiquement admissibles. Cette construction est entièrement parallèle et les champs $\underline{u}_D$, $\underline{u}_N$, $\lambda_N$ sont naturellement disponibles au sein des algorithmes. Le seul surcoût réside dans la reconstruction du champ statiquement admissible, mais ce coût est réduit puisque la reconstruction est parallèle. L'estimateur proposé dans [4] permet d'obtenir une borne garantie de la forme suivante :

$$ \|\underline{u}_{ex} - \underline{u}_D\|_{\mathbb{H}^{-1}, \Omega} \le \sqrt{\sum_s e_{CR_{\Omega}^{(s)}}^2 (\underline{u}_D, \underline{\hat{\sigma}}_N)} \quad (2) $$

Dans [6], une nouvelle borne garantie séparant ces deux sources d'erreur a été démontrée :

$$ \|\underline{u}_{ex} - \underline{u}_N\|_{\mathbb{H}^{-1}, \Omega} \le \sqrt{\sum_s e_{CR_{\Omega}^{(s)}}^2 (\underline{u}_N, \underline{\hat{\sigma}}_N)} + \sqrt{(r, z)} \quad (3) $$

où

— $\sqrt{(r, z)}$ est une mesure du résidu du solveur itératif : il s'agit du produit scalaire entre le résidu $r$ et le résidu préconditionné $z$.
— $\sqrt{\sum_s e_{CR_\Omega(s)}^2 (\underline{u}_N, \underline{\hat{\sigma}}_N)}$ se calcule en parallèle et est piloté par l'erreur de discrétisation sur les sous-domaines.

On illustre le résultat sur le problème élastique linéaire suivant (voir figure 1(a)). La structure a été décomposée en 16 sous-domaines et discrétisée en triangles. Le problème a été résolu en utilisant l'algorithme FETI et les champs statiquement admissibles ont été reconstruits par la technique EET avec un raffinement de type p+2 pour les résolution des problèmes locaux. Les résultats fournis par l'estimateur (3) sont présentés à la figure 1(b). Cette séparation permet d'envisager un nouveau critère d'arrêt du solveur itératif, du type stop quand $\sqrt{\langle r, z \rangle} \le \frac{1}{10} \sqrt{\sum_s e_{CR_\Omega(s)}^2 (\underline{u}_N, \underline{\hat{\sigma}}_N)}$, ce qui, sur l'exemple, permettrait d'arrêter le solveur à la sixième itération.

<figure><img src="image_2.png" /><figcaption>(a) Chargement du problème de référence (bleu) et du problème(b) Estimation d'erreur globale : séparation des adjoint (orange) et décomposition en 16 sous-domaines sources</figcaption></figure>

<caption>FIGURE 1 – Structure pré-fissurée</caption>

Dans [7], l'estimation d'erreur en cadre sous-structuré a été étendue aux quantités d'intérêt. Ces quantités, davantage susceptibles d'intéresser les ingénieurs par rapport à une information énergétique globale, sont définies à l'aide d'un extracteur [3] considéré comme une fonctionnelle du déplacement notée $\tilde{l}$ et supposée linéaire. On définit classiquement le problème suivant, dit problème adjoint :

$$ \text{Trouver } (\tilde{u}_{ex}, \underline{\hat{\sigma}}_{ex}) \in CA^0(\Omega) \times \tilde{SA}(\Omega) \text{ tel que } e_{CR_{\Omega}}(\tilde{u}_{ex}, \underline{\hat{\sigma}}_{ex}) = 0 \quad (4) $$

où $\tilde{SA}(\omega) = \{ \tau \in (L^2(\omega))_{\text{sym}}^{d \times d} ; \forall \nu \in CA^0(\omega), \int_{\omega} \tau : \underline{\varepsilon} = \underline{\varepsilon}(\nu) d\omega = \tilde{l}(\nu) \}$. La solution de ce problème existe et est unique. Le problème adjoint associé est résolu sur le même maillage que le problème de référence. Ainsi, la résolution des deux problèmes peut se faire simultanément avec un algorithme de type Block-Krylov [8]. De la même manière que dans [4] il est possible de construire les champs admissibles pour les deux problèmes. Les quantités notées avec un tilde proviennent du problème adjoint. Une borne classique sur les quantités d'intérêt peut être réécrite de manière à séparer les sources d'erreur :

$$ |I_{ex} - \tilde{l}(\underline{u}_D) - I_{HH,2}| \le \frac{1}{2}(\sqrt{r^T z} + e_{CR_\Omega}(\underline{u}_N, \underline{\hat{\sigma}}_N))(\sqrt{\tilde{r}^T \tilde{z}} + e_{CR_\Omega}(\tilde{u}_N, \tilde{\underline{\hat{\sigma}}}_N)) \quad (5) $$

$$ \text{avec } I_{HH,2} = \frac{1}{2} \int_{\Omega} (\underline{\hat{\sigma}}_N - \mathbb{H}\underline{\varepsilon}(\underline{u}_D)) : \mathbb{H}^{-1} : (\underline{\hat{\sigma}}_N + \mathbb{H}\underline{\varepsilon}(\tilde{u}_D)) d\Omega $$

Comme précédemment, cette séparation permet de définir un critère d'arrêt du solveur en lien avec la discrétisation : le solveur est piloté par un objectif de précision sur la quantité d'intérêt.

Afin d'améliorer la borne, nous proposons un enrichissement de la discrétisation du ou des sous-domaines contribuant majoritairement à l'erreur de discrétisation. Ces sous-domaines sont identifiés à l'aide des

$$ \text{mesures suivantes : } \eta^{(s)} = \frac{e_{CR_{\Omega(s)}} (\underline{u}_N^{(s)}, \underline{\hat{\sigma}}_N^{(s)})}{\sqrt{\sum_s e_{CR_{\Omega(s)}}^2 (\underline{u}_N^{(s)}, \underline{\hat{\sigma}}_N^{(s)})}} \text{ et } \tilde{\eta}^{(s)} = \frac{e_{CR_{\Omega(s)}} (\tilde{u}_N^{(s)}, \tilde{\underline{\hat{\sigma}}}_N^{(s)})}{\sqrt{\sum_s e_{CR_{\Omega(s)}}^2 (\tilde{u}_N^{(s)}, \tilde{\underline{\hat{\sigma}}}_N^{(s)})}} $$

Nous proposons un remaillage basé sur un raffinement hiérarchique qui induit des discrétisations imbriquées. Bien que créant une incompatibilité de maillage à l'interface, nous montrons qu'il est toujours possible de reconstruire de manière parallèle les champs requis sans conséquence sur la qualité
de l'estimation. Ce raffinement a été effectué sur le sous-domaine contributaire majoritaire qui est le 6 (qui porte la singularité du problème de référence et le chargement du problème adjoint). Nous obtenons les répartitions d'erreur présentées figure 2, où l'on constate que le niveau global d'erreur a diminué et que le sous-domaine 6 contribue autant que le sous-domaine 8 après raffinement. Enfin, la discrétisation de l'interface restant inchangée via ce processus, nous pouvons alors mettre en place des techniques de recyclage des espaces de Krylov où les directions de recherche générées lors d'une première résolution sont réutilisées lors de la suivante comme contrainte additionnelle accélérant ainsi la convergence du solveur.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Distribution d'erreur par sous-domaine</figcaption></figure>

Références

[1] P. Gosselet, C. Rey. *Non-overlapping domain decomposition methods in structural mechanics*, Archives of Computational Methods in Engineering, Vol 13, 515–572, 2006.

[2] P. Ladevèze, J.P. Pelle. *Mastering calculations on linear and non-linear mechanics*, Springer, 2006.

[3] M. Paraschivoiu, J. Peraire, A.T. Patera, *A posteriori finite element bounds for linear-functional outputs of elliptic partial differential equations*, Computer Methods in Applied Mechanics and Engineering, Vol. 150, 289–312

[4] A. Parret-Freaud, C. Rey, P. Gosselet, and F. Feyel. *Fast estimation of discretization error for FE problems solved by domain decomposition*, Computer Methods in Applied Mechanics and Engineering, Vol. 199, 3315–3323.

[5] F. Pled, L. Chamoin ; P. Ladevèze, *On the techniques for constructing admissible stress fields in model verification : performances on engineering examples*. Int Jal for Numerical Methods in Engineering 88 (5), 409–441, 2011

[6] V. Rey, C. Rey, P. Gosselet, *A strict error bound with separated contributions of the discretization and of the iterative solver in non-overlapping domain decomposition methods*, Computer Methods in Applied Mechanics and Engineering. Vol. 270, 293–303, 2013.

[7] V. Rey, P. Gosselet, C. Rey, *Strict bounding of quantities of interest in computations based on domain decomposition*, submitted.

[8] Y. Saad, *Analysis of augmented Krylov subspace methods*, Society for Industrial and Applied Mathematics, Vol 18, 435–449, 1997.

[9] C. Farhat, F.X. Roux, *Implicit parallel processing in structural mechanics*, Computational Mechanics Advances, Vol 2, 1–124, 1994.

[10] P. Le Tallec. *Domain-decomposition methods in computation mechanics*, Computational Mechanics Advances, Vol 1, 121–220, 1994.