Scalable parallel 3-D TEM inversion via rational approximation of the matrix exponential

Résumé technique : Inversion 3-D TEM parallèle et évolutive par approximation rationnelle de l'exponentielle de matrice

Énoncé du problème
L'électromagnétisme transitoire (TEM) est une technique cruciale en géophysique appliquée pour l'exploration minière, l'hydrogéologie et les études environnementales. Cependant, l'inversion de modèles de conductivité 3-D dans des environnements non triviaux est extrêmement exigeante sur le plan computationnel. Les approches traditionnelles reposent sur des schémas d'intégration temporelle implicites (par exemple, Euler implicite) pour résoudre les équations de Maxwell dans le domaine temporel. Ces méthodes organisent le flux de travail computationnel autour d'une marche temporelle, où le coût de la modélisation directe et de l'inversion est étroitement couplé au nombre de pas de temps requis pour résoudre la réponse transitoire. Alors que les gains de performance sur un seul cœur plafonnent en raison des limites de la loi de Moore, ces algorithmes intrinsèquement séquentiels peinent à évoluer efficacement sur le matériel parallèle moderne.

Méthodologie
Les auteurs proposent un cadre d'inversion novateur qui découple le coût computationnel du nombre de temps d'observation en reformulant le problème direct à l'aide de l'action d'une exponentielle de matrice.

1. Approximation rationnelle de l'exponentielle de matrice :
   Au lieu d'utiliser une marche temporelle, la solution transitoire $u(t)$ est exprimée comme l'action de l'exponentielle de matrice $\exp(-t M^{-1}K)$ sur un vecteur initial. Les auteurs emploient une approximation rationnelle quasi-optimale (RBA) de cette exponentielle, construite à l'aide de l'algorithme RKFIT. Cette approximation représente la solution comme une combinaison linéaire de systèmes linéaires décalés :
   $$u(t_j) \approx 2 \text{Re} \left( \sum_{i=1}^m \alpha_i^{(j)} (K - \xi_i M)^{-1} f \right)$$
   où $\{ \xi_i \}$ sont des pôles partagés entre tous les points temporels, et $\{ \alpha_i^{(j)} \}$ sont des coefficients précalculés. Crucialement, les matrices des systèmes décalés $(K - \xi_i M)$ sont indépendantes du temps.

2. Calcul parallèle de la solution directe et des sensibilités :
   La reformulation permet de résoudre les $m$ systèmes linéaires décalés de manière indépendante et parallèle. La dépendance temporelle est entièrement contenue dans les résidus (coefficients $\alpha$), ce qui signifie que le calcul des réponses directes et des sensibilités est effectivement indépendant du nombre de temps d'observation souhaités $K$.

   • Résolveur direct : Utilise des résolveurs directs parallèles (MUMPS) pour factoriser et résoudre les systèmes décalés.
   • Sensibilités (Jacobien) : Le Jacobien est calculé implicitement en utilisant le même ensemble de systèmes décalés. La dérivée de la matrice de masse est appliquée aux vecteurs de solution précalculés $g_i = (K - \xi_i M)^{-1}f$. Cela évite les récursions directes/adjointes séquentielles typiques des méthodes de marche temporelle.

3. Cadre d'inversion :

   • Optimisation : Le problème inverse est résolu à l'aide d'un cadre de Gauss-Newton. Le problème des moindres carrés linéarisé est résolu via l'algorithme LSQR, qui traite le Jacobien comme un opérateur implicite (sans matrice).
   • Régularisation : Les contraintes de régularité sont formulées à l'aide d'éléments de Raviart-Thomas (RT0) pour gérer des paramètres de modèle constants par morceaux sur des maillages non structurés, assurant une discrétisation cohérente du gradient.
   • Implémentation : Le code est écrit en Julia, utilisant le cadre Gridap.jl pour l'assemblage par éléments finis et une interface MPI personnalisée pour MUMPS afin de gérer les factorisations parallèles sur des architectures à mémoire partagée.

Contributions clés

• Changement de paradigme dans la conception algorithmique : L'article démontre un passage de l'adaptation de méthodes de marche temporelle séquentielles au matériel parallèle, à la conception de formulations mathématiques (RBA) qui exposent intrinsèquement la concurrence.
• Découplage des canaux temporels : Le coût computationnel est déterminé par le nombre de pôles $m$ (généralement faible, par exemple 16–21) plutôt que par le nombre de canaux temporels $K$ (souvent des centaines). Cela permet un échantillonnage temporel haute résolution sans augmenter la charge d'algèbre linéaire.
• Flux de travail d'inversion évolutive : En intégrant le parallélisme au niveau du modèle mathématique, les auteurs réalisent une inversion évolutive où le coût dominant (résolution des systèmes décalés) peut être distribué sur des rangs MPI.
• Publication en open source : L'implémentation est publiée sous forme de logiciel open source pour soutenir la recherche reproductible au sein de la communauté géophysique.

Résultats
La méthode a été validée sur un modèle synthétique 3-D d'environ 700 000 degrés de liberté, présentant des anomalies conductrices et résistives dans un demi-espace.

• Précision de reconstruction : L'inversion a permis de récupérer avec succès les emplacements et les formes générales des anomalies. Les blocs conducteurs ont été bien résolus, tandis que les blocs résistifs ont montré un contraste réduit, cohérent avec les limitations connues de la sensibilité TEM et de la régularisation par régularité.
• Performance :
  • Mise à l'échelle : Les expériences de mise à l'échelle forte ont démontré des performances quasi idéales jusqu'à 8 rangs MPI pour la factorisation et la substitution arrière. L'efficacité a légèrement diminué à 16 rangs en raison des effets NUMA (Non-Uniform Memory Access) sur la mémoire partagée et de la contention de la bande passante mémoire.
  • Mémoire : La consommation maximale de mémoire était d'environ 620 Go, nécessaire pour stocker simultanément les factorisations de 21 systèmes décalés.
  • Comparaison : Par rapport à la marche temporelle d'Euler implicite (qui nécessite 31 factorisations à valeurs réelles pour la même fenêtre temporelle), l'approche RBA nécessite 21 factorisations à valeurs complexes. Bien que l'arithmétique complexe soit plus coûteuse, la méthode RBA élimine le chemin critique de la marche temporelle séquentielle, offrant un potentiel supérieur de réduction du temps réel sur des architectures distribuées.

Signification et revendications
Les auteurs affirment que ce travail fournit une alternative évolutive à l'inversion TEM 3-D standard, adressant spécifiquement le goulot d'étranglement de la récursion temporelle. L'avantage principal n'est pas nécessairement une réduction du travail arithmétique sériel, mais la réduction du chemin critique parallèle et l'indépendance du coût vis-à-vis du nombre de canaux temporels.

L'article reconnaît les limitations actuelles, notant que le stockage simultané de multiples factorisations creuses crée une empreinte mémoire significative et que les architectures à mémoire partagée font face à des goulots d'étranglement de bande passante. Les auteurs suggèrent que les travaux futurs devraient se concentrer sur la transition vers des clusters de calcul haute performance (HPC) distribués et potentiellement le remplacement des résolveurs directs par des résolveurs itératifs préconditionnés pour améliorer davantage l'évolutivité et réduire l'utilisation de la mémoire. L'étude conclut que le cadre proposé récupère avec succès des structures 3-D et établit une base pour des inversions électromagnétiques à grande échelle plus efficaces.