Une Approche Multi-échelles pour des Calculs de Structures sur Ordinateurs à Architecture Parallèle

Préconditionneur de Neumann pour le problème primal

   Dans le cas de figure où tous les sous-domaines possèdent une partie de frontière à déplacement imposé (∂1Ω E 6= 0/), les condensées locales de Schur [SE ] sont toutes inversibles. La matrice du complément de Schur [S] étant l’assemblage des [S E], un préconditionneur envisageable est celui qui correspond aux problèmes de Neumann locaux : [S¯] −1 est l’assemblage des [S E ] −1 [De Roeck et Le Tallec 91], [Le Tallec 94]. L’application de [S¯] −1 à un résidu conduit à résoudre des problèmes de NeumannDirichlet (puisque ∂1Ω E 6= 0/) sous-domaine par sous-domaine pour déterminer le champ de déplacement bord correspondant [q E Γ ]. Le champ résultant [qΓ] n’étant plus continu, il est nécessaire de procéder à une opération de pondération (ou de prise de moyenne) des déplacements sur l’interface Γ issus de deux sous-domaines voisins (voire plus dans le cas de nœuds « coins » en 2D et 3D, et de nœuds « arêtes » en 3D). Pour traiter la phase de préconditionnement, il va donc être nécessaire d’avoir la factorisée de [K E ], matrice de rigidité du sous-domaine ΩE . Or, dans la méthode primale, seule la partie [K EE] est factorisée ; en particulier pour des problèmes de renumérotation, on sera conduit à stocker et factoriser deux matrices distinctes, ce qui peut être très pénalisant sur l’encombrement. De plus, la phase de préconditionnement conduit à pratiquer une montée-descente supplémentaire ; le coût d’une itération comme celui de l’initialisation est donc doublé en première approximation.

Décomposition en sous-structures et interfaces

  Comme précédemment, afin de rompre la globalité de la recherche d’un champ de déplacement et d’un champ de contrainte admissibles sur toute la structure Ω, on introduit une partition du milieu, mais cette fois-ci en deux entités distinctes : les sousstructures et les interfaces, [Ladevèze 87], [Ladevèze 96]. La vision du problème proposée ici consiste à considérer ces deux entités comme des structures à part entière, chacune possédant ses propres variables et les équations traduisant son propre comportement. Lorsque la structure étudiée résulte d’un assemblage de sous-structures, la partition s’introduit naturellement en utilisant les surfaces de jonction entre les diverses sous structures. Lorsqu’elle est composée d’une ou d’un nombre réduit de sous-structures massives, la partition peut se faire de façon artificielle au sein des sous-structures, par exemple, à partir de la seule géométrie, en utilisant une décomposition automatique [Farhat 88], [Malone 88], [Padovan et Kwang 91], [Al-Nasra et Nguyen 91], [Dagum 93], [Farhat et Lesoinne 93]. Cette partition peut alors répondre à des impératifs d’optimisation (réduction des communications entre processeurs, vitesse de convergence de l’algorithme) ainsi qu’à l’équilibrage des charges entre processeurs parmi lesquels seront distribuées les sous-structures.

Implantation — Synchronisme

  L’implantation de l’algorithme précédent a été réalisée dans le cadre du code éléments finis de type industriel CASTEM 2000, développé au CEA ([Verpeaux et al. 88]). Il s’agit d’un code orienté objets intégrant un pré et post-processeur, qui opère sur une base de donnée d’objets par une suite d’opérations élémentaires écrites dans le macrolangage de commande « Gibiane ». Les opérateurs élémentaires sont quant à eux développés dans le langage de programmation « Esope », [Verpeaux 93], transformé en Fortran par un pré-compilateur. Outre la gestion dynamique de la mémoire, l’environnement de programmation utilisé, « Gemat », permet un débordement programmé sur disque. Afin de tester l’algorithme LATIN sur une architecture parallèle de type MIMD à mémoire partagée ou distribuée, le code prototype spécialisé, un micro-CASTEM 2000,n’utilise que certaines des fonctionnalités d’origine (gestion dynamique des objets et débordement sur disque) et ne possède pas de langage de commande : il ne réalise que la partie résolution du problème, le pré-traitement (réalisation du maillage, décomposition en sous-structures et interfaces) et post-traitement étant réalisés sur une machine « hôte » qui supporte l’intégralité du code, [Millard 92]. Les données du calcul issus du pré-traitement ainsi que les résultats à post-traiter sont transférés par fichier pour pouvoir réaliser le traitement en « batch » sur des machines de production. La distribution des sous-structures et des interfaces parmi les processeurs de la machine se fait suivant la technique proposée dans [Lorong 94] : un groupe de sous structures et ses interfaces connectées sont attribués à chaque processeur. Certaines interfaces s’en trouvent dupliquées mais les calculs locaux sur ces dernières étant peu volumineux, cette technique n’est pas pénalisante. Les échanges d’information entre processeurs sont réalisés par le code, exécuté sur tous les processeurs, avec une technique d’échange de messages. Le code prototype utilise donc, suivant le cas, soit des bibliothèques natives (sur nCUBE 2 et iPSC 860),soit un « standard » d’échange de messages. Le choix de ce dernier s’est porté sur le système PVM (Parallel Virtual Machine), [Geist et al. 94], [Champaney et al. 96], mais les fonctionnalités requises sont similaires à celles offertes par toutes ces bibliothèques (MPI, NX-lib, . . .), figure 2.7.

Méthodes multi-grilles

   Historiquement initiées dans [Southwell 35] et [Fedorenko 64], ces méthodes se sont principalement développées pour les problèmes de mécanique des fluides à partir de 1970. Parmi les applications plus récentes au calcul de structures, citons par exemple [Parsons et Hall 90a], [Lebon et al. 89]. Leur principe ([Hackbush et Trottenberg 81]) repose sur le fait que les méthodes itératives classiques (gradient conjugué, Gauss-Seidel, Jacobi, relaxation, . . .) sont plus efficaces pour réduire les résidus à faible longueur de variation au cours des itérations, que pour les grandes longueurs de variation. Cela se traduit par une rapide décroissance de l’erreur lors des premières itérations puis une stagnation de celle-ci, ou du moins un taux de convergence diminué. Par exemple, considérons le problème d’une poutre en traction uniformément chargée, discrétisée en 31 éléments identiques (31 segments à 2 nœuds pour 31 degrés de liberté) dont on perturbe la solution de référence éléments finis par tous les modes de la matrice de rigidité (avec une même énergie de déformation sur chaque mode) comme l’illustre la figure 4.1.

Éléments finis hiérarchiques

   On souhaite tout d’abord déterminer une solution grossière sur un maillage grossier afin de trouver les zones à raffiner localement. Le but des éléments finis hiérarchiques est de permettre le calcul de la nouvelle solution sur le nouveau maillage tout en conservant le maximum d’informations issues du travail précédent, à savoir : une solution approchée grossière, et une matrice de rigidité grossière assemblée (et éventuellement factorisée) [Adjedj et Aubry 89]. Le principe de ces éléments finis hiérarchiques a été donné dans [Zienkiewicz et al. 83] : il consiste à remplacer la base éléments finis classique par une base dite hiérarchique. Le problème Kq = f posé sur la grille fine peut alors être réécrit en séparant les fonctions de base en deux groupes : celles associées aux éléments grossier (d’indice C) et les fonctions de base complémentaires (d’indice F).

Représentation du problème sur deux échelles

   Une représentation sur deux échelles qui peuvent être issues de deux modélisations différentes du même problème de départ, va nécessairement faire intervenir les quantités cinématiques (déplacements, déformations, . . .) et leurs quantités statiques duales (efforts, contraintes, . . .) sur chacune des deux modélisations. Bien évidemment, ces deux modélisations correspondent à deux échelles, dans le sens où les phénomènes qu’elles permettent d’appréhender dans la réponse de la structure possèdent des longueurs de variation différentes. Il faut alors permettre à ces deux modèles de « dialoguer » entre eux. Cela peut être réalisé par l’intermédiaire d’une paire d’opérateurs : suivant la terminologie des approches multi-grilles, [Hackbush et Trottenberg 81], [Le D

Conclusion

   Dans ce travail, nous nous sommes plaçés dans le cadre des techniques de décomposition de domaine, en vue de réaliser des calculs de structures sur ordinateurs à architecture parallèle. L’algorithme de base utilisé est issu des approches à grand incrément de temps lorsqu’elles sont couplées à une technique de sous-structuration. Ici, seul le comportement élastique linéaire des sous structures a été envisagé : il s’agit donc d’un cas de figure dégénéré, pour lequel seule la configuration finale nous intéresse, par rapport aux problèmes non-linéaires d’évolution pour lesquels la méthode LATIN a été initialement développée. À cause de cela, la dualité effort déplacement est substituée à celle qui est naturelle — entre effort et vitesse — dans les situations où le temps joue un rôle majeur (dynamique, problèmes d’évolution, etc.). Pour améliorer les performances de cette approche, une extension multi-échelles de cette dernière est mise en place. Elle consiste à appréhender la solution du problème de structure sur deux échelles distinctes. La première, l’échelle globale, est liée à une longueur caractéristique de la structure ; la deuxième, l’échelle fine, concerne les effets de longueur de variation beaucoup plus courte. En cherchant la solution successivement dans les deux échelles, et en conservant la résolution par éléments finis classique du problème à grande échelle, l’information se propage parmi toutes les sous-structures à chaque itération de l’algorithme. L’implantation dans le cadre du code éléments finis semi-industriel CASTEM 2000 a permis de montrer la faisabilité de l’approche. La réalisation de ce code prototype spécialisé a conduit à utiliser des calculateurs MIMD à mémoire distribuée avec un nombre de processeurs de l’ordre de 64. Comme le comportement du matériau est ici élastique linéaire, l’objectif est donc la réalisation de calculs de structures implicites de grande taille sur ces architectures parallèles. Dans ce type de situations, l’approche proposée est comparée à différentes démarches existantes de décomposition de domaine. L’influence de la taille des problèmes et des coûts relatifs des problèmes à grande échelle est illustrée sur des exemples à géométrie simple. Un aspect numérique important est la mise en œuvre d’une discrétisation traitant à égalité d’importance déplacements et efforts aux interfaces. En particulier, l’influence du choix de cette discrétisation sur la solution obtenue à convergence a été étudiée. Enfin, un aspect à souligner concerne le contrôle des paramètres de la méthode, c’est à dire la valeur de la longueur caractéristique des directions de recherche de l’algorithme itératif. Avec une approche mono-échelle, il est montré que sa valeur optimale est liée à la dimension de la structure ; avec l’utilisation de l’extension multiéchelles, elle devient cette fois-ci un paramètre d’interface, lié à la sous-structuration qui a été choisie. Dans la réalisation qui en a été faite ici, le code prototype traite essentiellement les cas de figure où les maillages éléments finis des deux échelles sont emboîtés, ce qui contraint la géométrie et donne peu de souplesse à l’utilisateur. Le traitement à partir de deux maillages indépendants serait donc un important gain en flexibilité. Une autre perspective à court terme est la possibilité de traitement d’interfaces non parfaites (comme par exemple le contact unilatéral qui a été présenté) à l’intérieur de l’approche multi-échelles. Bien entendu, l’objectif à plus long terme est le traitement de problèmes non linéaires tels que ceux engendrés par un comportement plastique ou viscoplastique du matériau, ou par une modélisation plus fine du comportement comme celle nécessitée, par exemple, par l’analyse des structures composites. C’est ce type de  problèmes qui entre dans le cadre du calcul intensif de structures, et qui nécessitera, de plus, le développement de techniques de contrôle des algorithmes et de la qualité des solutionsobtenues.

Table des matières

Introduction
1 Méthodes de décomposition de domaine 
1.1 Problème de référence et notations
1.2 Approches en déplacement
1.3 Approches duales
1.4 Préconditionnements
1.5 Bilan
2 Une approche à grand incrément de temps 
2.1 Décomposition en sous-structures et interfaces
2.2 Algorithme de résolution
2.3 Convergence
2.4 Exemples
3 Discrétisation
3.1 Cas d’une formulation en déplacement
3.2 Discrétisations et propriétés à convergence
3.3 Unicité
3.4 Implantation numérique
3.5 Exemples
4 Élaboration d’une approche multi-échelles 
4.1 Techniques à plusieurs échelles
4.2 Représentation du problème sur deux échelles
4.3 Couplage avec une décomposition en sous-structures
4.4 Algorithme itératif
4.5 Implantation numérique
5 Exemples, performances et comparaisons 
5.1 Influence des différentes échelles
5.2 Performances de l’approche proposée
5.3 Comparaisons avec la méthode FETI
Conclusion
Bibliographie
Annexes 1
A Intégration bord approchée, choix des points d’intégration
B Relation locale de comportement du problème 1D
C Algorithme du gradient conjugué
D Technique de lagrangien augmenté

Télécharger le rapport complet