# Décomposition de domaine pour le calcul des maçonneries

Alexandre Martin, Vincent Visseq, Paul Taforel, Frédéric Dubois

► To cite this version:

Alexandre Martin, Vincent Visseq, Paul Taforel, Frédéric Dubois. Décomposition de domaine pour le calcul des maçonneries. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722339

HAL Id: hal-01722339

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

Submitted on 3 Mar 2018

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
# Décomposition de domaine pour le calcul des maçonneries

Alexandre MARTIN <fnref n="1" />, Vincent VISSEQ <fnref n="2" />, Paul TAFOREL <fnref n="2" />, Frédéric DUBOIS <fnref n="2" />,<fnref n="3" />

**Résumé** — Une stratégie de décomposition de domaine appliquée à la méthode NSCD est testée pour le calcul d'ouvrages maçonnés. Précédemment introduite et validée pour l'analyse des milieux granulaires denses, la méthode Non Smooth Contact Domain Decomposition (NSCDD) est ici adaptée au contexte de la tenue mécanique des bâtiments. Une analyse de performance de cette méthode est mise en œuvre et comparée aux performances obtenues par parallélisation multitâches pour machines à mémoire partagée.

**Mots clés** — calcul parallèle, maçonneries, éléments discrets, dynamique non régulière, décomposition de domaine

## 1 Introduction

Les temps de calcul des simulations avec la méthode NSCD (Non Smooth Contact Dynamics [4]) sont généralement très longs en comparaison de calculs par éléments finis linéaires et même non linéaires. Le calcul parallèle offre une solution pour réduire le temps de restitution d'une simulation. Deux approches de parallélisation sont envisagées dans le logiciel LMGC90 [3] :

- le multi-threading, s'appuyant sur l'interface de programmation OpenMP pour tirer partie de la structure multi-cœurs des processeurs actuels. Elle permet à un programme de s'exécuter sur plusieurs cœurs de calcul. La contrainte majeure est que le programme ne peut s'exécuter que sur un nœud de calcul à la fois (mémoire partagée). Cette approche est disponible dans LMGC90 depuis plusieurs années [6, 7],

- le multi-processing, s'appuyant sur la librairie de passage de message MPI pour tirer partie de plusieurs ordinateurs connectés sur un réseau rapide. Le principe de la méthode consiste à exécuter plusieurs instances d'un même programme (processus) dialoguant entre elles. Cette méthode nécessite des modifications, parfois profondes, du code de calcul considéré. Les communications entre processus en sont l'un des facteurs limitant. Une méthode de décomposition de domaines (DDM), développée actuellement au LMGC [1, 8], a été adaptée à la simulation de structures maçonnées.

## 2 Méthode NSCD et parallélisation

### 2.1 Dynamique des contacts

#### 2.1.1 Equation de la dynamique

La dynamique de l'ensemble des corps considérés rigides va s'écrire classiquement :

$$ M \dot{\mathbf{V}} = \mathbf{F}_{\text{int}} + \mathbf{F}_{\text{ext}} + \mathbf{R} \quad (1) $$

avec $M = \bigcup_j M^j$ la matrice assemblée pour l'ensemble des corps $j$ des matrices $M^j = \begin{pmatrix} m^j I & 0 \\ 0 & J^j \end{pmatrix}$,
$\mathbf{V} = \bigcup_j \{\mathbf{v}^j, \omega^j\}^T$, $\mathbf{F}_{\text{ext}} = \bigcup_j \{\mathbf{f}^j_{\text{ext}}, \mathbf{m}^j_{\text{ext}}\}^T$ et $\mathbf{F}_{\text{int}} = \bigcup_j \{0, -\omega^j \times J^j \omega^j\}^T$. On constate que le système 1 est non-linéaire à cause du terme $\omega^j \times J^j \omega^j$. De par l'occurrence de chocs, ou

<footnote n="1">¹ LaMSID - Laboratoire de Mécanique des Structures Industrielles Durables, alexandre-externe.martin@edf.fr</footnote>

<footnote n="2">² LMGC, Laboratoire de Mécanique et Génie Civil, {vincent.visseq,paul.taforel,frederic.dubois}@univ-montp2.fr</footnote>

<footnote n="3">³ Laboratoire MIST, IRSN - CNRS - université de Montpellier 2</footnote>
d'autres événements non-réguliers, l'accélération ou la force de contact devront être comprises
au sens des mesures différentielles.

Par la suite on manipulera la forme intégrée sur un pas de temps, qui s'apparente à un bilan
de quantité de mouvement sur un pas de temps :

$$
\mathbb{M}(\mathbf{V}_{i+1} - \mathbf{V}_i) = \int_{t_i}^{t_{i+1}} \mathbf{F}_{\text{int}} + \mathbf{F}_{\text{ext}} dt + \mathbf{I} \quad (2)
$$

soit encore :

$$
\mathbf{V}_{i+1} = \underbrace{\mathbf{V}_i + \mathbb{M}^{-1} \int_{t_i}^{t_{i+1}} \mathbf{F}_{\text{int}} + \mathbf{F}_{\text{ext}} dt}_{\mathbf{V}_{\text{free}}} + \mathbb{M}^{-1} \mathbf{I} \quad (3)
$$

Le terme de vitesse libre **V**<sub>free</sub> est obtenu en intégrant la résultante et le moment des forces extérieures et le terme de rotation en utilisant un θ-schéma. **I** est l'impulsion des efforts de contact qui contient aussi bien des termes réguliers, que non réguliers.

2.1.2 Passage local-global

Le passage entre niveaux global et local, nécessaire à la résolution de ce type de problème, s'effectue à l'aide d'opérateurs linéaires $\mathbb{H}(q)$ et $^T\mathbb{H}(q)$ construits grâce à des considérations cinématiques et de dualité, et qui permettent de relier l'ensemble des vitesses relatives au contact $\mathcal{V}$ et l'ensemble des vitesses des corps $\mathbf{V}$ ou encore de calculer la résultante des impulsions de contact $\mathbf{I}$ à partir des impulsions de contact $\mathcal{I}$ :

$$
\mathbb{H}(q) : (\mathcal{I}^{\alpha}) \rightarrow \mathbf{I} = \sum_{\alpha} \mathbb{H}^{\alpha}(q) \mathcal{I}^{\alpha} \text{ et } {}^{T}\mathbb{H}(q) : \mathbf{V} \rightarrow (\mathcal{V}^{\alpha}) = \sum_{\alpha} {}^{T}\mathbb{H}^{\alpha}(q) \mathbf{V}
$$

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

Grace aux opérateurs de passage on peut réécrire l'équation de la dynamique sous une forme
condensée au contact :

$$
\mathbf{V}_{i+1} = \mathbf{V}_{\text{free}} + \mathbb{W}\mathcal{I} \quad (4)
$$

où on a introduit l'opérateur de Delassus $\mathbb{W} = {}^T\mathbb{H}(\mathbf{q}_m)\mathbb{M}^{-1}\mathbb{H}(\mathbf{q}_m)$.

2.1.3 Loi de contact frottant

On utilise classiquement une loi de contact en vitesse :

$$
\begin{cases}
\text{si } g > 0, \mathcal{I}_n = 0, \mathcal{V}_n \text{ quelconque} \\
\text{si } g \le 0, \mathcal{V}_n \ge 0, \mathcal{I}_n \ge 0 \text{ et } \mathcal{V}_n \cdot \mathcal{I}_n = 0
\end{cases}
$$

et une loi de frottement lorsque $\mathcal{I}_n \ge 0$ :

$$
\| \mathcal{I}_t \| \le \mu \mathcal{I}_n \left\{ \begin{array}{ll} \| \mathcal{I}_t \| < \mu \cdot \mathcal{I}_n \Rightarrow \mathcal{V}_t = 0 & \\ \| \mathcal{I}_t \| = \mu \cdot \mathcal{I}_n \Rightarrow \mathcal{V}_t = -\lambda \cdot \frac{\mathcal{I}_t}{\| \mathcal{I}_t \|} & \end{array} \right.
$$

où $\lambda > 0$ et $\mu$ est le coefficient de frottement à l'interface entre les deux corps.
### 2.1.4 Résolution

Si on considère les équations que doivent vérifier un contact $\alpha$ :

$$ egin{cases} 
\mathbf{V}_{i+1}^{\alpha} = \mathbf{V}_{free}^{\alpha} + \sum_{\alpha \neq \beta} \mathbf{W}^{\alpha\beta} \mathbf{I}^{\beta} + \mathbf{W}^{\alpha\alpha} \mathbf{V}^{\alpha} \\ 
\text{loi de contact } (\mathbf{I}^{\alpha}, \mathbf{V}^{\alpha}) = 0 
\end{cases} 
\quad (5) $$

On est capable de résoudre ce jeu d'équations si on considère les $\mathbf{I}^{\beta}$ connus. C'est le fondement de la méthode de Gauss-Seidel utilisée ici. On va résoudre successivement tous les contacts et itérer ces résolutions jusqu'à convergence. Ce qu'on résume ci dessous :

(0) Evaluating all the matrices $W^{\alpha\beta}$

$$ egin{cases} 
k = k + 1 \text{ (NLGS iteration)} \\ 
\alpha = \alpha + 1 \text{ (Contact index)} \\ 
\text{(a) Evaluating the right-hand side} \\ 
\mathbf{V}_{loc}^{\alpha} = \mathbf{V}_{free}^{\alpha} + \sum_{\beta<\alpha} \mathbf{W}^{\alpha\beta} \mathbf{I}^{k+1,\beta} + \sum_{\beta>\alpha} \mathbf{W}^{\alpha\beta} \mathbf{I}^{k,\beta} \\ 
\text{(b) Solving the local problem} \\ 
\text{Convergence test for } k = 0 \dots k_{max} 
\end{cases} $$

## 2.2 Parallélisation OpenMP

Une première solution pour permettre une gestion parallèle du code consiste à utiliser des directives OpenMP autour des boucles de calcul. Ainsi les boucles seront réalisées sur plusieurs "threads" simultanément. Cette technique ne nécessite pas de développement lourd et est assez paramétrable. Elle fonctionne comme la technique séquentielle si les opérations effectuées dans la boucle sont découplées.

De cette façon on peut paralléliser la boucle de contact. Même si les opérations ne sont pas strictement découplées cette technique fonctionne correctement [6].

(0) Evaluating all the matrices $W_{\alpha\beta}$

$$ egin{cases} 
k = k + 1 \text{ (NLGS iteration)} \\ 
\text{!$OMP PARALLEL PRIVATE (...) SHARED (...) ...} \\ 
\text{!$OMP DO ...} \\ 
\left[ 
\begin{array}{l} 
\alpha = \alpha + 1 \text{ (Contact index)} \\ 
\text{(a) Evaluating the right-hand side} \\ 
\mathbf{V}_{loc}^{\alpha} = \mathbf{V}_{free}^{\alpha} + \sum_{\beta<\alpha} \mathbf{W}^{\alpha\beta} \mathbf{I}^{k+1,\beta} + \sum_{\beta>\alpha} \mathbf{W}^{\alpha\beta} \mathbf{I}^{k,\beta} \\ 
\text{(b) Solving the local problem} \\ 
\text{!$OMP ENDDO} \\ 
\text{!$OMP END PARALLEL} \\ 
\text{Convergence test for } k = 0 \dots k_{max} 
\end{array} 
\right. 
\end{cases} $$

## 2.3 Parallélisation via la méthode NSCDD

La méthode de décomposition de domaine mise en œuvre ici (NonSmooth Contact Domain Decomposition [8]) consiste à sous-structurer le domaine et à :

- résoudre un problème de dynamique des contacts par sous-domaine, via l'algorithme itératif de Gauss-Seidel,

- introduire des conditions de recollement de la solution aux frontières des différents sous-domaines, permettant (à convergence) d'assurer la continuité des vitesses des corps d'interface, dans l'esprit de la méthode FETI.

Les principales différences des méthodes de décomposition de domaine adaptées à la simulation par éléments discrets, par rapport à celles relatives aux simulations par éléments finis, sont aux moins de trois types : le partitionnement en sous-domaines est réalisé à partir d'un graphe (évolutif) dont les sommets sont les corps et les segments sont les interactions et l'étape
<figure><img src="image_3.png" /><figcaption>Fig. 1 – Principales étapes de la méthode NSCDD pour un pas de temps.</figcaption></figure>

de décomposition doit être réalisée à chaque pas de temps (ou à une fréquence donnée), afin de prendre en compte les migrations des corps d'un sous-domaine à l'autre (cf. Fig. 1); la collection de corps traitée peut être fortement poly-disperse, pour reprendre une terminologie issue de la mécanique des milieux granulaires; l'algorithme de résolution de la dynamique des contacts est de type itératif et considère des relations non régulières entre inconnues de contacts.

## 3 Extension de la méthode au calcul de structures maçonnées

La méthode initialement développée pour la modélisation d'échantillons poly-disperses de disques (2D) ou de sphères (3D) a été étendue aux échantillons composés d'objets polyédriques afin de pouvoir être appliquée au calcul des structures maçonnées [2]. Une différence notable entre les maçonneries et les milieux granulaires provient de la faible migration d'objets entre sous domaines (point positif) et de la proportion des grains d'interface ayant une forte multipli-cité (point négatif).

Les structures maçonnées 3D peuvent en effet être considérées comme des assemblages de structures 2D; ce qui explique pourquoi la proportion des corps d'interface de ces structures est relativement faible en comparaison de celle rencontrée pour les milieux granulaires 3D (Ta-bleau 1).

<table><caption>Tableau 1 – Proportion des corps et des contacts d'interface pour des milieux granulaires et des structures maçonnées – cas 2D et 3D.</caption><thead><tr><td>Description du test</td><td>% particules d'interface</td><td>% contacts d'interface</td></tr></thead><tbody><tr><td>granulaire 2D (biaxial) 12 × 10³ disques</td><td>4.4%</td><td>10.7%</td></tr><tr><td>maçonnerie 3D 13 × 10³ polyèdres</td><td>4.6%</td><td>10.2%</td></tr><tr><td>granulaire 3D (triaxial) 50 × 10³ sphères</td><td>7.9%</td><td>17.1%</td></tr><tr><td>granulaire 3D (triaxial) 90 × 10³ sphères</td><td>9.1%</td><td>15.1%</td></tr></tbody></table>

Par ailleurs, la présence de corps de grande taille - utilisés pour la modélisation des dia-phragmes - par rapport à la dimension moyenne des blocs conduit à ce que la multiplicité as-sociée à ces corps soit très élevée en comparaison de celle que l'on trouve pour des échantillons granulaires classiques.

Enfin, il est également important de préciser que la phase de répartition des corps dans les sous-domaines peut être réalisée beaucoup moins fréquemment dans le cas des maçonneries que des milieux granulaires, la configuration de contact dans ce type de structure étant stable.

L'extension de la méthode NSCDD aux maçonneries a nécessité des adaptations dans la mise en œuvre numérique de la méthode de détection grossière du contact [5]. La détection grossière
du contact se déroule en deux étapes :

1. construction d'une liste de paires de contacteurs voisins par une méthode des boîtes utilisant la sphère englobante des contacteurs. Dans le cas des milieux granulaires, le réseau de contact occupe tout le volume de l'échantillon. Le domaine occupé par l'échantillon est découpé en boîtes qui *a priori* contiennent toutes une liste non vide de contacteurs. Dans le cas des maçonneries, le domaine englobant l'échantillon est globalement creux. Les boîtes situées à la périphérie du domaine concentrent la quasi-totalité des contacteurs, puisqu'elles intersectent l'enveloppe du bâtiment. Une nouvelle implémentation de la méthode des boîtes basée sur un stockage *sparse*, a donc été proposée. Les listes de contacteurs contenus dans les boîtes non vides sont stockées dans des listes chaînées plutôt que dans un tableau surdimensionné.

2. élimination des paires de contacteurs si il n'y a pas d'intersection des "axis aligned bounding box" (AABB) ni de *shadow-overlap*. Dans un premier temps on vérifie que les "ombres" sur les axes des deux contacteurs potentiellement en contact s'intersectent. Si cela n'est pas le cas, on élimine le contact de la liste. Par la suite, on applique une méthode de *shadow-overlap* qui repose sur la recherche d'un vecteur tel que la projection des "ombres" des deux contacteurs potentiellement en contact sur celui-ci ne s'intersectent pas. Dans le cas des maçonneries, on cherche ce vecteur parmi les normales aux faces.

Dans le cas du contact entre sphères, et pour une distance d'alerte bien choisie, la liste de paires de contacteurs sélectionnée fournit quasiment le réseau de contact et un faible nombre de paires de contacteurs sera éliminée lors de l'étape de détection fine. Dans le cas du contact entre polyèdres, l'anisotropie des contacteurs et la présence de contacteurs d'une taille bien plus élevée que la moyenne (typiquement les corps utilisés pour la modélisation des diaphragmes) imposent d'utiliser un filtrage supplémentaire.

Dans le cadre de la NSCDD, la mise à jour de la boîte englobant chaque contacteur polyèdrique est réalisée en parallèle sur chaque processeur, lors de la mise à jour de la configuration de détection. Les boîtes englobantes sont ensuite envoyées sur un processeur hôte afin d'effectuer la détection grossière.

# 4 Résultats numériques

L'analyse des performances, en termes de parallélisation de l'algorithme NSCD, est quantifiée par la suite à l'aide des quantités usuelles : le speed-up $S_p$, qui est le rapport du temps de calcul de la version mono-processeur au temps de calcul de la version parallèle sur $N_p$ processeurs, et le temps de restitution, qui est le temps de calcul perçu par le modélisateur. Afin d'illustrer le comportement des deux approches de parallélisation, deux études sont effectuées sur une structure maçonnée de petite taille à deux étages, représentée Fig. 2.

<figure><img src="image_4.png" /><figcaption>Fig. 2 – Structure maçonnée R+1 utilisée pour les études numériques.</figcaption></figure>

Pour les calculs menés en décomposition de domaine, il est nécessaire, en plus du nombre
de processus, de définir une grille de partitionnement ($n_x, n_y, n_z$). La figure 3 présente les partitionnements étudiés par la suite.

<figure><img src="image_5.png" /><figcaption>Fig. 3 – Partitionnements proposés pour l'étude de la méthode NSCDD. Grilles de partitionnement ($n_x, n_y, n_z$); 2 sous-domaines : (2, 1, 1); 4 sous-domaines : (2, 2, 1); 8 sous-domaines : (2, 2, 2).</figcaption></figure>

## 4.1 Mise sous poids propre d'une structure R+1

Dans cette première étude, l'influence du nombre de processus (1, 2, 4 et 8) sur le temps de restitution (Fig. 4) ainsi que sur le speed-up est analysé pour un chargement de type mise sous poids propre (activation de la gravité).

Malgré un temps séquentiel plus faible de la version OpenMP, avec compilateur propriétaire Intel (Ifort), vis-à-vis de la version DDM-MPI, avec compilateur open source (Gfortran), les temps de restitution pour 8 processeurs sont équivalents. Les speed-up des deux méthodes restent confondus jusqu'à 4 processus, puis divergent sensiblement pour 8 processus.

<figure><img src="image_6.png" /><figcaption>Fig. 4 – Temps de restitution en fonction du nombre de processus (a); speed-up en fonction du nombre de processus (b). Mise en évidence du coût de la gestion du parallélisme (overhead), pour les approches OpenMP et MPI, sur les temps de restitution (Figure 4(a)) et calcul du speed-up prenant en compte une référence mono-domaine avec ou sans gestion du parallélisme (Figure 4(b)).</figcaption></figure>

## 4.2 Tassement différentiel d'une structure R+1

Dans cette deuxième étude, l'influence du nombre de processus de la méthode NSCDD (égal au nombre de sous-domaines), sur le speed-up, est analysé pour un chargement de type tassement différentiel, définit Fig. 5(a) et Fig. 5(b).

Le nombre d'itérations par pas de temps pour les différentes sous-structurations (Fig. 6(a)) suit les mêmes tendances générales. Cependant, le nombre d'itérations total est sensiblement supérieur pour les structures décomposées en sous-domaines (de l'ordre de 20%).

Les interpénétrations maximale (Fig. 7(a)) et moyenne (Fig. 7(b)) présentent des différences modérées entre les versions mono-domaine et multi-domaines, ces valeurs d'interpénétrations restant très faibles devant les dimensions des blocs de la structure maçonnée.
<figure><img src="image_7.png" /><figcaption>Fig. 5 – Géométrie de la structure R+1 avec (en bleu foncé) la fondation pilotée en vitesse (a); profil de vitesse imposée en fonction du temps (b) : les deux zones gris foncé correspondent à la mise sous poids propre et à une vitesse de descente constante du plancher rigide, la zone gris clair représente un chargement de type tassement différentiel régularisé.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6 – Nombre d'itérations par pas de temps pour les différentes découpes en sous-domaines (1, 2, 4 et 8) (a); speed-up en fonction du nombre de sous-domaines (b).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7 – Évolution des interpénétrations maximale (a) et moyenne (b) au cours de la simulation pour les différentes découpes en sous-domaines.</figcaption></figure>

Le motif de fissuration en fin de calcul, sensiblement identique pour toutes les simulations réalisées, est présenté Fig. 8. L'obtention d'un même profil de fissuration pour l'ensemble des partitionnements, couplé à un speed-up proche du speed-up linéaire, permet de conclure à une bonne adaptation de la méthode de décomposition de domaine au calcul des maçonneries.
<figure><img src="image_10.png" /><figcaption>Fig. 8 – Motifs de fissurations de la structure maçonnée en fin de tassement différentiel et vitesses des éléments discrets suivant z.</figcaption></figure>

# 5 Conclusion

Les résultats de l'étude de performance montrent que la méthode de décomposition de domaine est adaptable aux simulations de structures maçonnées. Des différences entre les deux approches de parallélisation proposées sont mises en évidence. Ce type de stratégie doit permettre à l'avenir la modélisation de structures plus imposantes tout en maîtrisant les durées de restitution des simulations.

Remerciements :

Ce travail a pu être réalisé grâce à OSEO-FEDER (projet Degrip).

Travail effectué avec le support du Centre de compétences en calcul haute performance de la région Languedoc-Roussillon HPC@LR, équipé d'un calculateur hybride IBM, financé par la région Languedoc Roussillon, l'Europe et l'Université Montpellier 2 Sciences et Techniques.

# Références

[1] Alart, P., Iceta, D. and Dureisseix, *A nonlinear domain decomposition formulation with application to granular dynamics*, Computer Methods in Applied Mechanics and Engineering, 59-67, 2012.

[2] F. Dubois, A. Martin, *Rapport scientifique - projet DEGRIP*, 2012.

[3] F. Dubois, M. Jean, M. Renouf, R. Mozul, A. Martin et M. Bagnéris, *LMGC90*, 10ème Colloque National en Calcul des Structures, Giens, 2011.

[4] M. Jean, *The non-smooth contact dynamic method*, Computer methods in applied mechanics and engineering, 177:235-257, 1999.

[5] F. Radjaï, F. Dubois, *Modélisation numérique discrète des matériaux granulaires*, Hermès - Lavoisier, 2010.

[6] M. Renouf, P. Alart, *Solveurs parallèles pour la simulation de systèmes multi-contacts*, Revue européenne des éléments finis, vol.13, 691-702, 2004.

[7] M. Renouf, F. Dubois, P. Alart, *Numerical investigations of fault propagation and forced-fold using a non smooth discrete element method*, Revue européenne de mécanique numérique, vol.15, 549-570, 2006.

[8] V. Visseq, A. Martin, D. Dureisseix, F. Dubois, P. Alart, *Distributed Nonsmooth Contact Domain Decomposition (NSCDD) : algorithmic structure and scalability*, 21th Domain Decomposition International Conference, Rennes - France, 2012.