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

Riepilogo Tecnico: Inversione 3-D TEM Parallela Scalabile tramite Approssimazione Razionale dell'Esponenziale di Matrice

Enunciato del Problema
L'Elettromagnetismo Transitorio (TEM) è una tecnica critica nella geofisica applicata per l'esplorazione mineraria, l'idrogeologia e gli studi ambientali. Tuttavia, l'inversione di modelli di conduttività 3-D in ambienti non banali è computazionalmente onerosa. Gli approcci tradizionali si basano su schemi di integrazione temporale implicita (ad esempio, Eulero all'indietro) per risolvere le equazioni di Maxwell nel dominio del tempo. Questi metodi organizzano il flusso di lavoro computazionale attorno alla marcia temporale, dove il costo della modellazione diretta e dell'inversione è strettamente accoppiato al numero di passi temporali necessari per risolvere la risposta transitoria. Poiché i miglioramenti delle prestazioni single-core hanno raggiunto un plateau a causa delle limitazioni della Legge di Moore, questi algoritmi intrinsecamente sequenziali faticano a scalare in modo efficiente sull'hardware parallelo moderno.

Metodologia
Gli autori propongono un nuovo framework di inversione che disaccoppia il costo computazionale dal numero di tempi di osservazione riformulando il problema diretto utilizzando l'azione di un esponenziale di matrice.

1. Approssimazione Razionale dell'Esponenziale di Matrice:
   Invece dell'integrazione temporale, la soluzione transitoria $u(t)$ è espressa come l'azione dell'esponenziale di matrice $\exp(-t M^{-1}K)$ su un vettore iniziale. Gli autori impiegano un'approssimazione razionale quasi-ottimale (RBA) di questo esponenziale, costruita utilizzando l'algoritmo RKFIT. Tale approssimazione rappresenta la soluzione come una combinazione lineare di sistemi lineari traslati:
   $$u(t_j) \approx 2 \text{Re} \left( \sum_{i=1}^m \alpha_i^{(j)} (K - \xi_i M)^{-1} f \right)$$
   dove $\{ \xi_i \}$ sono poli condivisi tra tutti i punti temporali e $\{ \alpha_i^{(j)} \}$ sono coefficienti precalcolati. Crucialmente, le matrici dei sistemi traslati $(K - \xi_i M)$ sono indipendenti dal tempo.

2. Calcolo Parallelo Diretto e delle Sensibilità:
   La riformulazione consente di risolvere i $m$ sistemi lineari traslati in modo indipendente e parallelo. La dipendenza temporale è interamente contenuta nei residui (coefficienti $\alpha$), il che significa che il calcolo delle risposte dirette e delle sensibilità è effettivamente indipendente dal numero di tempi di osservazione desiderati $K$.

   • Solutore Diretto: Utilizza solutori diretti paralleli (MUMPS) per fattorizzare e risolvere i sistemi traslati.
   • Sensibilità (Jacobian): Il Jacobiano è calcolato implicitamente utilizzando lo stesso insieme di sistemi traslati. La derivata della matrice di massa viene applicata ai vettori di soluzione precalcolati $g_i = (K - \xi_i M)^{-1}f$. Questo evita le ricorsioni dirette/adiointe sequenziali tipiche dei metodi di integrazione temporale.

3. Framework di Inversione:

   • Ottimizzazione: Il problema inverso è risolto utilizzando un framework Gauss-Newton. Il problema dei minimi quadrati linearizzato è risolto tramite l'algoritmo LSQR, che gestisce il Jacobiano come operatore implicito (senza matrice).
   • Regolarizzazione: I vincoli di regolarità sono formulati utilizzando elementi Raviart–Thomas (RT0) per gestire parametri del modello costanti a tratti su mesh non strutturate, garantendo una discretizzazione coerente del gradiente.
   • Implementazione: Il codice è scritto in Julia, utilizzando il framework Gridap.jl per l'assemblaggio agli elementi finiti e un'interfaccia MPI personalizzata per MUMPS per gestire le fattorizzazioni parallele su architetture a memoria condivisa.

Contributi Chiave

• Cambiamento di Paradigma nella Progettazione Algoritmica: Il documento dimostra un passaggio dall'adattamento di metodi di marcia temporale sequenziali all'hardware parallelo, alla progettazione di formulazioni matematiche (RBA) che espongono intrinsecamente la concorrenza.
• Disaccoppiamento dei Canali Temporali: Il costo computazionale è determinato dal numero di poli $m$ (tipicamente piccolo, ad esempio 16–21) piuttosto che dal numero di canali temporali $K$ (spesso centinaia). Ciò consente un campionamento temporale ad alta risoluzione senza aumentare il carico di algebra lineare.
• Workflow di Inversione Scalabile: Incorporando il parallelismo a livello del modello matematico, gli autori raggiungono un'inversione scalabile in cui il costo dominante (risoluzione dei sistemi traslati) può essere distribuito tra i rank MPI.
• Rilascio Open-Source: L'implementazione è rilasciata come software open-source per supportare la ricerca riproducibile nella comunità geofisica.

Risultati
Il metodo è stato validato su un modello sintetico 3-D con circa 700.000 gradi di libertà, caratterizzato da anomalie conduttive e resistive in un mezzo semi-infinito.

• Accuratezza di Ricostruzione: L'inversione ha recuperato con successo le posizioni e le forme generali delle anomalie. I blocchi conduttivi sono stati ben risolti, mentre i blocchi resistivi hanno mostrato un contrasto ridotto, coerente con i noti limiti della sensibilità TEM e della regolarizzazione per regolarità.
• Prestazioni:
  • Scalabilità: Gli esperimenti di scaling forte hanno dimostrato prestazioni quasi ideali fino a 8 rank MPI sia per la fattorizzazione che per la sostituzione all'indietro. L'efficienza è degradata leggermente a 16 rank a causa degli effetti NUMA (Non-Uniform Memory Access) a memoria condivisa e della contesa della larghezza di banda di memoria.
  • Memoria: Il consumo di memoria di picco è stato di circa 620 GB, necessario per memorizzare le fattorizzazioni di 21 sistemi traslati simultaneamente.
  • Confronto: Rispetto all'integrazione temporale Eulero all'indietro (che richiede 31 fattorizzazioni a valori reali per la stessa finestra temporale), l'approccio RBA richiede 21 fattorizzazioni a valori complessi. Sebbene l'aritmetica complessa sia più costosa, il metodo RBA elimina il percorso critico sequenziale della marcia temporale, offrendo un potenziale superiore per la riduzione del tempo reale su architetture distribuite.

Significato e Affermazioni
Gli autori affermano che questo lavoro fornisce un'alternativa scalabile alla standard inversione TEM 3-D, affrontando specificamente il collo di bottiglia della ricorsione temporale. Il vantaggio principale non è necessariamente una riduzione del lavoro aritmetico seriale, ma la riduzione del percorso critico parallelo e l'indipendenza del costo dal numero di canali temporali.

Il documento riconosce le limitazioni attuali, notando che la memorizzazione simultanea di multiple fattorizzazioni sparse crea un significativo footprint di memoria e che le architetture a memoria condivisa affrontano colli di bottiglia nella larghezza di banda. Gli autori suggeriscono che il lavoro futuro dovrebbe concentrarsi sulla transizione verso cluster di calcolo ad alte prestazioni (HPC) distribuiti e potenzialmente sulla sostituzione dei solutori diretti con solutori iterativi precondizionati per migliorare ulteriormente la scalabilità e ridurre l'uso di memoria. Lo studio conclude che il framework proposto recupera con successo le strutture 3-D e stabilisce una base per inversioni elettromagnetiche più efficienti e su larga scala.