Une bibliothèque Python pour le développement de stratégies de calcul avancées

Programmation orientée objet pour les éléments finis et les méthodes de résolution associées

La programmation orientée objet offre plus de modularité, de ré-utilisabilité et de flexi-bilité qu’une programmation procédurale. De plus en plus d’outils de simulation numérique, académiques ou industriels, reposent sur le paradigme objet, notamment en langage C++. Les concepts de base de la programmation orientée objet peuvent se résumer à l’encap-sulation, l’héritage et le polymorphisme [Meyer, 1988]. Cependant de nombreux patrons de conception peuvent être constitués à partir de ces simples concepts [Wolfgang, 1994 ; Gamma et al., 1999]. Ces patrons de conception (ou Design Patterns) répondent à des problématiques récurrentes en conception logiciel. [Heng et Mackie, 2009] ont proposé des Design Patterns pour la conception d’un code éléments-finis. Nous proposons une brève revue des conceptions choisies dans la littérature.
Dans le cadre de la simulation numérique par la méthode des éléments-finis, des re-cherches ont donc été menées sur la façon de concevoir un code modulaire qui réponde aux forts besoins d’évolutivité de la méthode (formulation de nouveaux éléments, nouvelles équations constitutives, nouvelles méthodes de résolution, …). Des classes de base qui dé-crivent les entités d’un modèle éléments finis (degré de liberté, nœud, élément, maillage) ont été introduites par [Miller, 1988]. La formulation du problème en est indépendante et peut former une autre classe [Peskin et Russo, 1988]. Dans un cadre non-linéaire matériau, le point de Gauss et le matériau peuvent avantageusement faire partie des objets de base [Me-nétrey et Zimmermann, 1993]. [Dubois-Pèlerin et Pegon, 1998] se sont concentrés sur la conception des objets supervisant le modèle éléments finis pour l’analyse non-linéaire. Une classe Problem contient le modèle, une méthode de continuation et une méthode de résolu-tion. La méthode de résolution peut faire appel à des classes dérivées d’une classe abstraite (i.e. dont certaines méthodes ne sont pas implémentées) NonLinearSolver ou d’une classe abstraite LinearSolver selon les caractéristiques du problème. Il est alors aisé d’ajouter des classes NewtonRaphson, ModifiedNewtonRaphson ou BFGS qui implémentent la classe NonLinearSolver. De la même façon l’ajout de solveurs linéaires est simplifié par les classes abstraites qui définissent les entrées/sorties et les méthodes des solveurs linéaires et non-linéaires.
La problématique d’amélioration de la performance numérique a été abordée par [Pantale, 2005] à travers la parallélisation de l’assemblage dans un code orienté objet, ayant sensible-ment la même structure de classe que celle présentée plus tôt. L’encapsulation des données dans les classes de manière hiérarchisée (du degré de liberté au maillage complet) donne en effet une indépendance entre les éléments qui permet de paralléliser très directement le calcul des opérateurs élémentaires, via un protocole de multi-threading à mémoire partagée comme OpenMP [Dagum et Menon, 1998]. La parallélisation de la résolution peut se faire égale-ment par l’implémentation de méthodes de décomposition de domaine. Dans une bibliothèque existante (logiciel ZéBuLoN [Northwest Numerics, 2013]), [Gosselet, 2003] a ainsi proposé l’ajout d’une classe PARALLEL_FORMULATION. Cette classe dérive d’une classe IT_S_OPERATOR qui réalise les opérations classiques d’une boucle d’un solveur itératif ITER_SOLVER (produit matrice/vecteur, produit scalaire, préconditionnement, projection). La formulation parallèle repose sur une description des sous-domaines et de leurs interfaces par des objets spécifiques. Le traitement de chaque sous-domaine est affecté à un processeur, ce qui assure la localisation des données. Les degrés de liberté d’interface sont également stockés au niveau local dans une classe BOUNDARY. Une représentation des interfaces locales entre les domaines deux à deux (par connectivité) est privilégiée pour la communication entre les sous-domaines voisins. Les opérateurs d’interface (complément de Schur, préconditionneur, scaling, second membre condensé) sont aussi représentés par une classe qui dérive de la classe B_OPERATOR. Grâce à une instance d’une classe EXCHANGER qui propose des méthodes de communication (envoie et réception de message), la résolution itérative du problème d’interface (par un algorithme de type Krylov) peut être effectuée en parallèle. Les objets “compléments de Schur” disposent ensuite d’une méthode pour la localisation de la solution à l’intérieur des sous-domaines. Cette implémentation orientée objet de la méthode de décomposition de domaine autorise la définition de plusieurs types d’interface (primale, duale ou mixte) et l’ajout de différents opérateurs comme les préconditionneurs, les scalings, les problèmes grossiers.
En ce qui concerne la réduction de modèle, [Aquino, 2007] a proposé un implémenta-tion orientée objet de la POD. Une classe POD gère le sous-espace (snapshots, matrice de corrélation, base réduite POD). Le modèle réduit est également représenté par une classe ReducedModel qui réalise l’intégration des opérateurs élémentaires dans le sous-espace réduit plutôt que de projeter le système assemblé. La classe ReducedModel hérite, comme la classe FullFEModel du modèle complet, d’une classe BasePDE, ce qui rend possible le passage d’un modèle complet à un modèle réduit, et permet d’ajouter d’autres types de modèles réduits.
Une implémentation orientée objet de la méthode XFEM est proposée dans [Bordas et al., 2007]. D’une part, les classes de base d’un code éléments finis sont modifiées, ceci notamment pour intégrer l’ajout de degrés de liberté à certains nœuds. Le calcul de la matrice d’inter-polation et de la matrice des dérivées des fonctions d’interpolation des éléments est adapté pour prendre en compte l’enrichissement. Au niveau du problème global, l’enrichissement mène parfois à des dépendances linéaires dans la matrice de raideur qui sont supprimées par une méthode spécifique. D’autre part, des classes supplémentaires sont créées pour gérer les fonctions d’enrichissement (EnrichmentFunction) et la description des entités physiques et géométriques nécessitant un enrichissement tels que les fissures, les trous ou les interfaces entre matériaux. Une classe particulière pour les règles d’intégration numérique est proposée (classe IntegrationRule) ce qui permet de dériver plusieurs méthodes de la littérature [Be-lytschko et al., 2009].

Conception de la bibliothèque ICAFE

La conception de la bibliothèque ICAFE répond à un cahier des charges concernant le type de méthode qu’elle doit permettre de développer. Sans pouvoir tout implémenter à court terme, il est nécessaire d’anticiper les contraintes qui pourraient empêcher l’extension rapide vers les différents types de méthodes identifiés.

Les méthodes et leurs contraintes

L’étude bibliographique des stratégies de résolution avancées en Chapitre 2 a permis de mettre en lumière quelques approches actuelles mais il n’est bien sur pas envisageable de construire un code de recherche qui les intègre toutes. Au regard des travaux de cette thèse, des perspectives qui sont ouvertes, et des activités du laboratoire en termes de développement de méthodes numériques [Passieux et al., 2013b ; Gogu et al., 2010] certaines méthodes, et leur combinaison, sont privilégiées. Nous présentons ces méthodes et les contraintes associées en termes de fonctionnalités d’un code de recherche :
– Décomposition de domaine : différentes méthodes de décomposition de domaine sans recouvrement ont été présentées en Chapitre 2. Plusieurs types de raccords doivent pouvoir être définis afin d’étudier leur influence pour un problème donné, les précondi-tionneurs et les problèmes grossiers les plus adaptés. Il peut être également intéressant de définir une décomposition de domaine imbriquée (multi-niveaux) pour traiter effi-cacement des problèmes fortement multi-échelles.
– Réduction de modèle par projection : ces méthodes nécessitent de pouvoir construire une base réduite et de réaliser une projection de Galerkin du système d’équations. L’in-troduction de cette fonctionnalité, dans un solveur, ne présente pas de difficulté par-ticulière. L’intégration réduite des équations d’équilibre doit cependant être possible pour améliorer les performances de la réduction de modèle comme cela est envisagé pour les stratégies PBAMR et PBAMR-DD (voir sous-section 2.2.2 et Chapitre 4).
– Méthodes X/GFEM : l’enrichissement d’interpolation par les méthodes X/GFEM offre une grande flexibilité pour augmenter le degré d’interpolation ou introduire une description implicite de la géométrie. Pour cela, il faut prévoir l’ajout de degrés de li-berté aux nœuds et leur association avec des fonctions d’enrichissement. Les fonctions d’enrichissement doivent pouvoir faire appel à une description implicite de la géomé-trie (e.g. fonctions level set). Cette fonctionnalité permettrait d’approfondir l’étude de l’approche par enrichissement en post-flambement local présentée en Annexe A.
– Méthodes multi-grilles : dans ces méthodes, plusieurs discrétisations sont associées à la même zone géométrique, et celles-ci ne sont pas toujours imbriquées et ne coïn-cident pas forcément entre elles sur leurs contours. Il est donc important de pouvoir associer plusieurs maillages à une zone géométrique et de pouvoir construire les opéra-teurs de prolongation et de restriction qui permettent de transférer l’interpolation d’un champ d’un maillage à un autre. Certaines méthodes multi-grilles (Arlequin [Ben Dhia, 1998]) font évoluer un maillage raffiné local par dessus un maillage grossier (suivi d’un phénomène local, propagation de fissure) [Ben Dhia et Jamond, 2010 ; Passieux et al., 2013b].