# FastWaveED_GPU3D

GPU-accelerated 3D elastodynamic solver for anisotropic elastic media with full moment-tensor sources (earthquakes, explosions, collapses, etc.)
Implemented in CUDA C with MATLAB support for input, compilation, and postprocessing.
Single-GPU optimized.


**Version:** 1.3  
**Date:** 26 July 2025  
**Author:** Yury Alkhimenkov, Massachusetts Institute of Technology

![alt text](figs/lune_wavefield_diagram.png)


---

## Usage

### Input
Generated by MATLAB script `vizme_FastBiotED_GPU3D_OK_2024_final.m` (material parameters, moment tensor, etc.).

### Output
Seismograms in binary format (e.g., `*_Rec_y2.res`), read and visualized by MATLAB.

### Run steps
1. Open MATLAB and navigate to folder `FastWaveED_GPU3D_v1_DC_Earthquake_FULL_OK4_final`.
2. Execute:
   ```matlab
   vizme_FastBiotED_GPU3D_OK_2024_final
   ```
   This will:
   - Define model and parameters
   - Generate input files
   - Compile CUDA kernel, e.g.:
     ```bash
     nvcc -arch=sm_89 -O3 -DNBX=12 -DNBY=12 -DNZ=7 \
          -DOVERX=0 -DOVERY=0 -DOVERZ=0 \
          -DNPARS1=10 -DNPARS2=13 -DNPARS3=15 \
          FastWaveED_GPU3D_v1.cu
     ```
   - Run:
     ```
     Linux: bash filename
     Windows: a.exe
     ```
   - Visualize outputs in MATLAB

---

## Requirements

- Windows OS
- MATLAB (in system PATH)
- CUDA Toolkit ≥ 10.0
- NVIDIA GPU, compute capability ≥ 8.9

---
VERSION 1
Section 4.4 — An efficient source modeler and wave propagation solver (reformulated, drop‑in)

\subsection{An efficient source modeler and wave propagation solver}

Scope. A research‑grade toolchain for 
(i) moment‑tensor prototyping on the lune and 
(ii) 3‑D elastic wave propagation on GPUs. 
Ready to run for this paper’s use cases; fast; not a black‑box product.

Interfaces. 
A MATLAB driver orchestrates parameter sweeps, compiles the CUDA kernels with \texttt{nvcc}, launches runs, and collects seismograms/snapshots. 
Configuration via MATLAB structs and compile‑time macros; some changes require editing the CUDA sources.

Physics and numerics (as used here). 
Linear elastic, isotropic media; 
first‑order velocity–stress system on a staggered Cartesian grid; 
second‑order finite differences in space/time; 
Ricker source time function; 
point moment‑tensor injection; 
perfectly matched layers; 
single‑GPU execution; 
single or double precision. 
Governing equations and BCs as stated in §4.2.

Moment-tensor modeler. Moment tensors are generated on the lune; distances are measured by spherical arc length \citep{TapeTape2019}. 
We use explicit [11,12,13,22,23,33] [11,12,13,22,23,33] ordering when listing components to avoid ambiguity with Voigt. 
The modeler enumerates DC, ISO, CLVD, and mixed cases and emits tensors directly consumable by the GPU solver.

What is implemented but not relied on here. 
Hooks for 3‑D heterogeneity and anisotropy (orthorhombic stiffness in code) exist but are not exercised or validated in this paper; they are off by default. 
Multi‑GPU domain decomposition is not included in the release used here; any 8× figures shown elsewhere are ideal projections, not measured multi‑GPU runs.

Performance posture. 
Kernels are memory‑bandwidth bound and tile fields to coalesce global loads/stores; 
shared memory is used for stencil reuse; 
PML variables are stored contiguously. 
Effective throughput and runtimes are reported in §5 and Table~\ref{tab:results-setups}; 
reproducibility details (GPU models and precisions) are stated there.

I/O. Minimal, binary outputs for seismograms and optional 3‑D snapshots at user‑defined strides; 
no heavy post‑processing in the solver.

Limitations. 
Research code; sparse error handling; no packaged installer; 
changing stencil/PML thickness/field layout or adding a new source type may require CUDA edits and recompilation. 
No out‑of‑core execution; grid size is constrained by device memory.

Bottom line. 
The toolchain is fit for iterative MT studies and forward modeling at the scales demonstrated here, 
with measured speedups over CPU baselines. 
It remains research‑style software rather than a turnkey production system.

VERSION 2
The current implementation solves the coupled elastodynamic–poroelastic equations using a first-order finite-difference time-domain (FDTD) scheme on structured grids. 
The code targets NVIDIA GPUs via CUDA C and achieves high throughput. 
Both isotropic and orthorhombic stiffness tensors are supported. 
A flexible source injection scheme allows moment tensor sources (including full non-DC terms) and pressure sources. 
Material heterogeneity is fully supported.

Grid generation, model setup, and source specification are handled via MATLAB. 
A helper routine (mat2cuda) assembles input data and generates the appropriate kernel configuration. 
Kernel compilation is invoked from MATLAB and uses nvcc via system calls. 
Simulations are launched from the same interface. 
Output is stored in binary format for post-processing.

The solver supports absorbing boundaries via convolutional PML. 
Ricker wavelets are used for time functions. 
Temporal and spatial discretizations are controlled manually. 
No adaptive meshing or dynamic load balancing is implemented. 
GPU memory use is flat; full 3D grids are loaded to device memory. 
Maximum grid sizes depend on available memory (see Table X).

The code is performant and suitable for controlled experiments, parameter studies, and iterative inversion under user supervision. 
However, it is not a fully encapsulated production tool. 
Manual code edits are sometimes required to change model size, PML configuration, or physics flags. 
Grid sizes, material parameters, and kernel constants are currently hardcoded or passed at compile time.


## Citation

Cite both the software and the paper.

**Software (APA)**  
Alkhimenkov, Y. (2025). *FastWaveED_GPU3D* (Version 1.3) [Computer software]. Massachusetts Institute of Technology. GitHub. https://github.com/…

**Paper (APA example)**  
[All authors]. (2025). [Paper title]. *Geophysical Journal International*. https://doi.org/…


---

## License

Released under the [GNU General Public License v3.0](LICENSE).  
See the `LICENSE` file in this repository for full terms.