Utilisation de l'IBM Blue Gene/P de l'IDRIS :Le projet SEISCOPE

Modélisation numérique de la propagation des ondes sismiques en domaines temporel et fréquentiel : application à l?imagerie de l?intérieur de la Terre (SEISCOPE)

Responsable du projet : Stéphane Operto¹

Participants : H. Ben Hadj Ali¹, R. Brossier², V. Etienne¹, L. Giraud³, A. Haidar⁴, J. Virieux².

¹ Geoazur CNRS-IRD-UNSA-OCA, ² LGIT-UJF, ³ INRIA Toulouse, ⁴ CERFACS

Ce projet a pour objectif le développement de méthodes de modélisation 3D des ondes acoustiques et élastiques dédiées à des applications d?imagerie sismique par inversion du champ d?onde complet. La modélisation des ondes est effectuée dans les domaines temporel ou fréquentiel avec des méthodes aux différences finies et aux éléments finis de type Galerkin discontinu. Dans le domaine fréquentiel, une méthode de décomposition en domaines fondée sur une approche algébrique par complément de Schur et un solveur hybride direct/itératif est utilisée pour résoudre le système d?équations linéaires résultant de la discrétisation de l?équation d?onde acoustique 3D. Pour résoudre les équations de l?élastodynamique en 3D, une méthode éléments finis de type Galerkin discontinu sur maillage tétraédrique non structuré permet de tirer profit d?une adaptation locale de la taille des cellules et de l?ordre d?interpolation des solutions au sein de chaque cellule (adaptativité hp) pour traiter des milieux géologiques complexes impliquant de forts contrastes lithologiques. L?ensemble de ces méthodes de modélisation ont été ou sont en cours d?implémentation au sein de codes d?inversion. En 2009, une première application synthétique de la méthode d?inversion a été effectuée sur Babel où une technique de compression de données fondée sur l?encodage et l?assemblage des sources sismiques est utilisée pour réduire le coût des simulations multi-sources. Notre ambition en 2010 est d?appliquer cette approche à un jeu de données réelles enregistrées sur le champ pétrolier de Valhall.

L?imagerie du sous-sol a de nombreuses applications dans les domaines du génie civil (détection de cavités, instabilités gravitaires), de l?aléa sismique (étude des effets de site), de l?environnement (stockage de déchets), de l?énergie (prospection pétrolière et gazifère) et de l?étude fondamentale des phénomènes géodynamiques. Dans ce contexte applicatif, nous développons des méthodes d?imagerie sismique par inversion du champ d?onde complet (FWI : Full Waveform Inversion) dans le but d?extraire de manière aussi complète que possible l?information contenue dans les données sismiques. Les modèles du sous-sol sont reconstruits avec une résolution théorique maximale d?une demi-longueur d?onde et sont paramétrés par un ou plusieurs paramètres physiques gouvernant la propagation des ondes sismiques (vitesses de propagation des ondes de compression et de cisaillement, densité, atténuation, anisotropie). L?inversion est formulée en domaine fréquentiel, où des algorithmes d?imagerie multirésolution sont implémentés par inversions successives de fréquences croissantes. Un ingrédient essentiel des méthodes d?inversion est la modélisation de la propagation des ondes sismiques. Nous développons des méthodes de modélisation en domaines fréquentiel et temporel. Dans le domaine fréquentiel, nous utilisons des méthodes de décomposition en domaines fondées sur l?utilisation de solveurs direct et hybride pour résoudre le système d?équations linéaires résultant de la discrétisation de l?équation d?onde. En domaine temporel, V. Etienne a développé une méthode par éléments finis de type Galerkin discontinu sur maillage tétraédrique non structuré permettant une adaptation locale du maillage en termes de taille de cellules et d?ordre d?interpolation au sein de chaque cellule (adaptativité hp). L?ensemble de ces moteurs de modélisation ont été ou sont en cours d?implémentation au sein de codes d?inversion 3D pour la reconstruction de milieux 3D acoustiques et élastiques.

Code FWM3D_DDM_A : Modélisation en domaine fréquentiel fondée sur un solveur hybride direct/itératif (développeurs : H. Ben Hadj Ali, L. Giraud, A. Haidar, F. Sourbier).

Un code de modélisation 3D des ondes acoustiques, fondé sur une méthode de décomposition en domaines et une approche par complément de Schur, a été installé sur Blue Gene/P. Le domaine de calcul est décomposé en sous domaines sur grille Cartésienne et chaque sous-domaine est attribué à un processus MPI. Le système du complément de Schur fournissant les solutions aux noeuds d?interface est résolu en parallèle avec un solveur itératif GMRES tandis que le solveur direct MUMPS est utilisé pour la factorisation de matrices d?impédance locales assemblées sur chaque processeur. Une fois que les solutions aux noeuds d?interface ont été calculées avec GMRES, les solutions aux noeuds intérieurs aux domaines sont efficacement calculées par des produits matrice-vecteur et des phases de substitutions locales. Une efficacité supérieure ou égale à 1 est obtenue jusqu?à un nombre critique de sous-domaines au delà duquel la précision du préconditionneur du système de Schur est trop dégradée. Le code est écrit en Fortran 90 et en MPI.

Code FWM3D_DG_ELAST : Résolution des équations de l?élastodynamique en 3D par méthode éléments finis de type Galerkin discontinu temps-espace (développeur : V. Etienne).

Un code massivement parallèle de modélisation des ondes en milieu élastique 3D a été développé et installé sur Babel par V. Etienne, thésard à Geoazur. La discrétisation spatiale est fondée sur une méthode aux élements finis de type Galerkin discontinu tandis que la discrétisation temporelle est fondée sur un schéma classique aux différences finies d?ordre 2. Le temps de calcul est optimisé en tirant profit de l?adaptativité hp qui permet d?adapter localement la taille des cellules mais aussi l?ordre d?interpolation des solutions au sein de chaque cellule. Les conditions absorbantes sont des C-PML. Le domaine de calcul est représenté sur des maillages tétraédriques non structurés construits avec le mailleur Tetgen. Le partitionnement du maillage pour l?exécution parallèle du code est effectué avec METIS. Le code est écrit en Fortran 90 et repose sur une parallélisation MPI efficace due au caractère local des éléments finis permettant de limiter le volume de communication.

Code FWT3D_A (développeur : H. Ben Hadj Ali): Inversion du champ d?onde complet et encodage de sources

Un code d?inversion du champ d?onde complet utilisant une technique de compression de données par encodage de sources a été installé sur Babel par H. Ben Hadj Ali. La modélisation des ondes est calculée en domaine fréquentiel avec le solveur hybride direct/itératif implémenté dans FWM3D_DDM_A. La totalité des sources de l?acquisition sismique sont encodées et sommées avant la simulation pour se ramener à des problèmes de simulation mono-source.

Modélisation en domaine fréquentiel fondée sur un solveur hybride direct/itératif

Nous avons évalué la scalabilité et la précision du code FWM3D_DDM_A par comparaison avec une méthode classique aux différences finies temps-espace. Dans l?approche temporelle, le champ monochromatique est extrait de la simulation en domaine temporel en calculant une transformée de Fourier discrète dans la boucle en temps. Des problèmes impliquant 30 millions d?inconnues ont été traités sur 1024 processeurs en mode SMP. Nous avons testé un nouveau préconditionneur du système de complément de Schur qui a permis une diminution notable du nombre d?itérations de GMRES particulièrement en présence d?atténuation intrinsèque du milieu.

Modélisation dans des milieux élastiques 3D complexes.

La précision et la scalabilité du code FWM3D_DG_ELAST ont été validées sur différents exemples numériques de complexité croissante (milieu homogène infini, demi-espace). Une simulation sur le benchmark EUROSEISTEST représentatif d?un bassin sédimentaire situé dans la région de Thessalonique en Grèce a été effectuée et validée avec une solution calculée par la méthode des éléments spectraux (solution fournie par E. Chaljub (LGIT)) (Figure 1). Le maillage tétraédrique utilisé dans la méthode DG comprend 16,3 millions d'éléments pour un total de 1,31 milliards d?inconnues. La simulation a été calculée sur 4096 processeurs de Babel en mode VN. Le temps calcul était de 7 h.

Imagerie d?un modèle acoustique 3D par inversion des formes d?onde avec encodage de sources

H. Ben Hadj Ali a testé son code d?inversion 3D sur un jeu de données synthétiques représentatif d?une zone d?étude en prospection pétrolière (modèle SEG/EAGE overthrust). La taille du modèle est 20 km x 20 km x 5 km. Trois fréquences (3.5, 5 et 7 Hz) ont été successivement inversées en utilisant une technique d?encodage de sources (Figure 2). L?objectif en 2010 sera d?appliquer cette même technique à un jeu de données enregistré sur le champ pétrolier de Valhall.

Figure 1 :

Carte du mouvement maximum du sol calculée pour le benchmark EUROSEISTEST représentatif du bassin sédimentaire de Thessaloniki en Grèce. (à gauche) Solution calculée avec la méthode éléments finis de type Galerkin discontinu. (à droite) Solution calculée par la méthode des élements spectraux (carte fournie par E. Chaljub (LGIT)). La solution DG a été calculée sur 4096 processeurs en mode VN. Le nombre d?inconnues en DG est de 1,31 milliards. Le temps calcul est de 7 h. (Figure V. Etienne).

Figure 2 :

Imagerie du modèle overthrust par inversion du champ d?onde complet. A gauche: coupes verticale et horizontale du modèle exact. A droite: coupes reconstruites par inversion pour une acquisition réduite à une seule source assemblée. Trois fréquences (3.5, 5 et 7 Hz) ont été inversées en effectuant 200 itérations par fréquence. Notez le bruit résiduel associé aux interférences entre les sources encodées de la source assemblée (Figure H. Ben Hadj Ali).

• H. Ben Hadj Ali, S. Operto and J. Virieux Efficient 3D frequency-domain Full Waveform Inversion (FWI) with phase encoding. Extended abstracts. 71th EAGE Conference, Amsterdam 2009, P004.
• H. Ben Hadj Ali, S. Operto and J. Virieux Three-dimensional frequency-domain full waveform inversion with phase encoding. Extended abstracts. SEG Houston 2009 Annual meeting, p. 2288-2292.
• V. Etienne, J. Virieux, N. Glingsky and S. Operto Seismic modelling with Discontinuous Galerkin Finite-Element method - Application to Large Scale 3D Elastic Media. Extended abstracts. 71th EAGE Conference, Amsterdam 2009, P131.
• V. Etienne, J. Virieux and S. Operto A massively parallel time-domain Discontinuous Galerkin method for 3D elastic wave modeling. Extended abs. SEG Houston 2009 Annual meeting, p. 2657-2661.
• J. Virieux, S. Operto, H. Ben Hadj Ali, R. Brossier, V. Etienne, F. Sourbier, L. Giraud, A. Haidar, Seismic wave modeling for seismic imaging, The Leading Edge, 28, 538 - Special section: Seismic modeling, 2009.
• J. Virieux and S. Operto. An overview of full waveform inversion in exploration geophysics, Geophysics, 74(6), WCC127-WCC152, 2009.