Rapid Earthquake-to-Tsunami Waveform Generation via Large-Scale Multi-GPU FFT Convolution Applied to the Cascadia Subduction Zone

Technical Summary: Rapid Earthquake-to-Tsunami Waveform Generation via Large-Scale Multi-GPU FFT Convolution

Problem Statement
Data-driven early warning systems for earthquakes and tsunamis rely on large ensembles of rupture scenarios and their resulting waveforms to perform Bayesian inference and forecasting. However, generating these datasets via repeated high-fidelity numerical simulations of seismic and tsunami dynamics is computationally prohibitive. Conventional approaches require solving partial differential equations (PDEs) for tens of thousands of earthquake realizations, a process that is too slow for rapid ensemble evaluation. Specifically, for the Cascadia Subduction Zone (CSZ), existing workflows struggle to generate the necessary large-scale waveform ensembles within practical timeframes without sacrificing physical fidelity.

Methodology
The authors propose a pipeline that exploits the Linear Time-Invariant (LTI) structure inherent in both the elastodynamic equations (seismic wave propagation) and the acoustic-gravity wave equations (tsunami propagation). Instead of solving PDEs for every realization, the method precomputes impulse responses and reduces the source-to-waveform mapping to two consecutive convolution operators.

1. Precomputation of Impulse Responses:

   • Elastic Green's Functions: The first operator maps fault slip-rate functions to vertical seafloor displacement. Using the frequency-wavenumber (FK) method, the authors precompute Green's functions for 963 triangular subfaults. These functions are constructed offline and stored in the frequency domain.
   • Acoustic-Gravity Adjoint Responses: The second operator maps seafloor motion to offshore pressure and tsunami waveforms. Rather than performing forward solves for every seafloor grid point (which would be $N_m \approx 2.4$ million), the authors utilize an adjoint formulation. This allows them to compute the response by performing $N_q$ (64) adjoint solves, one for each observation location, effectively reversing the input-output relationship to reduce computational cost.

2. Block Toeplitz Formulation:
   The time-discretized dynamics result in two consecutive block lower-triangular Toeplitz matrix-vector multiplications (matvecs). The complete map transforms slip-rate functions ($s$) to seafloor motion ($m$) and then to observation waveforms ($q$):
   $$q = F_q D_t F_G s$$
   where $F_G$ and $F_q$ are the Toeplitz operators representing the elastic and acoustic-gravity dynamics, respectively, and $D_t$ is a temporal difference operator.

3. Distributed FFT-Accelerated Pipeline:
   To evaluate these massive convolutions efficiently, the authors implement a distributed algorithm using FFTMatvec.

   • Domain Decomposition: The large seafloor grid is partitioned across multiple GPUs. The elastic operator ($F_G$) is partitioned by output rows (seafloor points), while the adjoint operator ($F_q$) is partitioned by input columns (matching seafloor points).
   • In-Device Processing: The intermediate seafloor field is kept resident on the GPUs between the two convolution steps, avoiding the I/O bottleneck of writing and reading the full space-time field.
   • Frequency-Domain Acceleration: Local Toeplitz matvecs are evaluated using Fast Fourier Transforms (FFTs). The response kernels are transformed into the frequency domain once during setup and remain in GPU memory. The convolution is performed as element-wise multiplication in the frequency domain, followed by an inverse FFT.

Key Contributions

• Multi-GPU, FFT-Accelerated Source-to-Waveform Generation: The paper introduces a pipeline that formulates the earthquake-to-tsunami mapping as two consecutive block Toeplitz matvecs evaluated via distributed FFTs. This approach eliminates the need to solve wave equations for every realization and keeps intermediate fields on-device.
• Scalability to Full CSZ Scale: The authors demonstrate the pipeline on a full-scale Cascadia Subduction Zone configuration involving 963 subfaults, 2,416,530 seafloor grid points, 64 observation locations, and 256 timesteps.
• Memory Management: The implementation successfully manages the massive memory requirements of the frequency-domain kernels (approximately 9.45 TiB aggregate GPU memory) by distributing the data across 64 GPUs within a single NVL72 domain.

Results
The pipeline was evaluated on the TACC Vista system using NVIDIA GH200 and GB200 GPUs.

• Performance: On 64 GB200 GPUs, the pipeline generates waveforms for one rupture in 24 ms once the response operators are resident.
• Speedup: This represents a 21.8× speedup over a distributed PyTorch baseline implementation using direct time-domain convolution.
• Scaling:
  • Strong Scaling: Increasing from 8 to 64 GPUs yielded speedups of 6.5× (GH200) and 6.8× (GB200), with parallel efficiencies of 81.2% and 85.0%, respectively.
  • Weak Scaling: The system maintained high efficiency (95.1% on GH200 and 95.8% on GB200) as the problem size scaled from 8 to 64 GPUs.
• Ensemble Generation: At the measured rate, applying the operators to an ensemble of $10^4$ ruptures requires approximately 4 minutes, compared to 87 minutes for the PyTorch baseline.

Significance
The paper claims that this pipeline provides a practical route for generating large-scale waveform ensembles based on high-fidelity models. By reducing the source-to-waveform map to precomputed convolutions and leveraging distributed FFT acceleration, the method enables the evaluation of large rupture ensembles within minutes. This capability is critical for data-driven tsunami modeling and forecasting, particularly for Bayesian inverse problems where rapid generation of synthetic waveforms is required to constrain source parameters and predict tsunami heights. The work demonstrates that it is possible to bypass expensive-to-solve PDEs for every realization while maintaining the physical fidelity required for early warning systems.