DeepGPR

Official Website / Documentation: https://songc0a.github.io/DeepGPR/

DeepGPR provides a wave propagation module for PyTorch, designed for applications such as Ground Penetrating Radar (GPR) imaging and inversion. Its core concepts are derived from Deepwave. You can use it to perform both forward modeling and backpropagation—thereby enabling the simulation of wave propagation to generate synthetic data—as well as for Full Waveform Inversion (FWI). Furthermore, you can integrate this wave propagation functionality into a larger operational pipeline—incorporating various wavelets, loss functions, and other components—to achieve end-to-end forward and reverse propagation, powered by automatic differentiation and our high-performance operators.

Features

Supports 2D and 3D forward modeling of Maxwell's equations—via the Finite-Difference Time-Domain (FDTD) method—for both single and multiple excitation scenarios.

Gradients of the output receiver data can be computed with respect to model parameters (relative permittivity, conductivity), the initial wavefield, and source amplitudes.

Automatically extends the physical model into CPML, with an independent thickness for each boundary. Supply only air and the target region; material gradients retain the input model shape.

The compute backend can run on CUDA GPUs or on CPU. The CPU backend is implemented in C and is selected automatically when device='cpu'.

The FDTD spatial finite-difference order can be selected with fdtd_order=2, 4, or 8 (default: 2).

The FWI gradient mode can be selected with mode=2 or mode=3. mode=2 keeps the previous Ez-only gradient behavior, while mode=3 uses Ex, Ey, and Ez forward/adjoint electric-field contributions for relative permittivity and conductivity gradients.

Supports techniques such as checkpointing, DDP, and the utilization of CPU memory to minimize GPU memory consumption, thereby enabling the execution of large-scale models.

System Requirements

• OS: Linux, Windows, and macOS for CPU execution; Linux and Windows for CUDA execution
• Environment: Python 3.8+, CUDA Toolkit for CUDA execution
• Libraries: torch, numpy, scipy, matplotlib
• Hardware: NVIDIA GPU with sufficient VRAM for CUDA execution; CPU execution works without a GPU.

Start

Before CUDA use, you must ensure that you have an NVIDIA graphics card and have installed a CUDA-enabled version of PyTorch. For CPU use, install a CPU build of PyTorch and include a compiled deepgpr_cpu shared library in src/DeepGPR/lib.

DeepGPR can then be installed using

  pip install DeepGPR

A Small Forward Modeling Test

import torch
import DeepGPR
import matplotlib.pyplot as plt

# Set up the parameters and models
device=torch.device("cuda" if torch.cuda.is_available() else "cpu")
dx=0.02  # Or [dx, dy, dz], for example [0.02, 0.015, 0.01]
dt=3e-11
nt=2000
er = torch.ones(100, 100,1) * 2  
er[50:,:]=5
se = torch.zeros_like(er)  
er.requires_grad_()
source_location=torch.tensor([[[10,10,0]]],device=device,dtype=torch.int)
receiver_location=torch.tensor([[[10,90,0]]],device=device,dtype=torch.int)
freq=2e8
peak_time = 1 / freq
source_amplitudes = torch.zeros((1,nt,1),device=device)
source_amplitudes[0,:,0]=DeepGPR.wavelet.ricker(
    freq, nt, dt, peak_time, device=device
)


#forward modeling
r = DeepGPR.compute(
    device=device, dx=dx, dt=dt, 
    source_amplitudes=source_amplitudes,
    source_location=source_location, 
    receiver_location=receiver_location, 
    er=er, se=se,
    fdtd_order=2
)

(r[-1]**2).sum().backward()

_, ax = plt.subplots(1, 2, figsize=(10, 3))
ax[0].plot(r[-1].detach().flatten().cpu().numpy())
ax[0].set_title("Receiver data")
ax[1].imshow(er.grad.detach())
ax[1].set_title("Gradient")
plt.show()

There are more examples in the ./examples.

Source Wavelets

Wavelets are available from the DeepGPR.wavelet module. Every function returns a one-dimensional PyTorch tensor and accepts optional dtype and device arguments.

ricker = DeepGPR.wavelet.ricker(freq, nt, dt, peak_time, device=device)
gaussian = DeepGPR.wavelet.gaussian(freq, nt, dt, peak_time, device=device)
derivative = DeepGPR.wavelet.gaussian_derivative(
    freq, nt, dt, peak_time, device=device
)
morlet = DeepGPR.wavelet.morlet(
    freq, nt, dt, peak_time, cycles=3.0, device=device
)
burst = DeepGPR.wavelet.sine_burst(
    freq, nt, dt, peak_time, cycles=3.0, device=device
)

DeepGPR.ricker(...) remains available as a backward-compatible alias.

The following figures present representative 2D and 3D full-waveform inversion (FWI) examples. For each case, the true model, initial model, and inverted result are shown to evaluate the reconstruction performance of the proposed method.

2D FWI Result

The 2D example illustrates the inversion performance on a two-dimensional subsurface model. The comparison between the true model, initial model, and inverted result shows that the proposed method can effectively recover the main structural features from the initial model.

True Model	Initial Model	Inverted Result

3D FWI Result

The 3D example demonstrates the applicability of the proposed method to three-dimensional full-waveform inversion. For visualization, the figures below show the central slice of the 3D model, including the true model, the initial model, and the inverted result. The comparison indicates that the proposed method can reconstruct the dominant subsurface structures in the 3D case and improve the model consistency relative to the initial model.

Model Type	Central Slice of 3D Model
True Model
Initial Model
Inverted Result

compute Interface Documentation

compute is a core function for 3D/2D Finite-Difference Time-Domain (FDTD) forward modeling, primarily designed for Ground Penetrating Radar (GPR) and electromagnetic wave propagation. It fully supports backpropagation (e.g., for Full Waveform Inversion, FWI) utilizing PyTorch's autograd engine.

📝 Function Signature

def compute(device, dx=None, dt=None, 
            source_amplitudes=None,
            source_location=None, 
            receiver_location=None, 
            er=None, se=None, mr=None, 
            E=None, H=None, PML=None,
            pmlthick=10, source_direction=2, reciever_direction=2,
            model_gradient_sampling_interval=1,
            wavefield_storage_dtype=torch.float32,
            use_async_offload=False,
            fdtd_order=2,
            mode=2,
            debug=False,
            print_parameters=False,
            save_forward_wavefield_path=None):

Set print_parameters=True on an ordinary call to print the normalized configuration and memory estimate immediately before the native solver starts:

result = DeepGPR.compute(
    device="cuda:0",
    # Other model, source, and acquisition arguments...
    print_parameters=True,
)

📥 Input Parameters

1. Basic Physics & Grid Parameters

2. Medium Model Parameters

Supply only the physical model, including any air layer and the target region. Do not include PML in eps_r, sigma, or mu_r. For 2D simulations use (nx, ny) or (nx, ny, 1).

Every compute call replicates the current material values at each boundary outward by pmlthick. For [px0, px1, py0, py1, pz0, pz1], the internal grid is Nx = nx + px0 + px1, Ny = ny + py0 + py1, Nz = nz + pz0 + pz1. In 2D, Nz = 1. Sources and receivers use zero-based input model coordinates; the solver adds [px0, py0, pz0] internally without modifying your tensors. PML thickness can exceed the physical model size.

After loss.backward(), eps_r.grad and sigma.grad have exactly their respective input shapes, including air cells. Update these physical tensors with your external FWI optimizer; the next call regenerates PML from the updated model. If air must remain fixed, apply your physical air mask in the optimizer.

Migration: remove old manually padded PML cells and subtract the old low-face padding from acquisition indices. Keep actual air layers. Existing states generated on a different computational grid must be regenerated.

Dimension Key: nx, ny, and nz represent the number of grid cells along the X, Y, and Z axes, respectively.

3. Source & Receiver Setup

This section defines the geometric observation system (coordinates) and the excitation waveforms.

Core Shape Definitions:

• nstep: Number of shots/batches (independent simulation tasks running in parallel).
• nsr: Number of sources per single simulation.
• nrx: Number of receivers per single simulation.
• nt: Total number of time steps to simulate.

4. Boundary Conditions & Optimization

4.1 FWI Gradient Mode

mode only changes how the model gradients are accumulated during backpropagation:

• mode=2 (default): Saves Ez in E_saved and computes relative permittivity/conductivity gradients from Ez only.
• mode=3: Saves Ex, Ey, and Ez in E_saved and computes relative permittivity/conductivity gradients from all three electric-field components. This is intended for complete 3D Maxwell FWI. The receiver component is not changed by this option.

When eps_r or sigma requires gradients, mode=2 is restricted to 2D Ez-TM modeling; use mode=3 for 3D gradients. Source-waveform gradients are supported. Gradients with respect to mu_r are not currently implemented and are rejected explicitly.

4.2 Discrete Adjoint Gradient

The backward solver applies the exact reverse-mode transpose of each executed operation in reverse order: receiver sampling, source injection, electric CPML, electric update, magnetic CPML, and magnetic update. The derivative transpose is applied to the material-weighted field cotangent, so heterogeneous update coefficients and anisotropic grid spacing are handled by the executed discrete operator. Every electric and magnetic CPML auxiliary state has a separate cotangent recurrence on all six faces.

CPML is treated as a fixed numerical boundary during backward. Boundary coefficient averages are detached, and gradients are cropped to the physical model; replicated PML sensitivities are not summed into model edges. All physical cells remain eligible for material gradients, including the first cell at a low face. PML materials and coefficients are rebuilt on the next forward call. Thus a finite-difference check must keep boundary values fixed (or freeze the extended PML and its coefficients explicitly); perturbing model edges and regenerating PML also changes a numerical boundary that this FWI gradient intentionally holds fixed.

The material-gradient formulation in DeepGPR was informed in part by the differentiable FDTD implementation in TIDE, particularly its treatment of the discrete Maxwell electric-field update in gradient computation. We gratefully acknowledge the TIDE project and its authors for this work.

Use model_gradient_sampling_interval=1 and wavefield_storage_dtype=torch.float32 for a directional derivative check in the physical model region. Temporal subsampling and lower-precision storage deliberately approximate the gradient. Run the gradient-check notebook after rebuilding the native ABI 6 libraries with deepgpr_supports_external_pml. See tests/test_external_pml.py for physical edge and checkpoint regression checks.

4.3 GPU-native block INT8 history

Only the sampled forward state used by the material adjoint is compressed. The live Ex/Ey/Ez/Hx/Hy/Hz and CPML states remain float32. The executed discrete update requires E^n and R^n, where E^(n+1) = ca E^n + cb R^n; consequently mode=2 stores compressed Ez/Rz and mode=3 stores compressed Ex/Ey/Ez and all three corresponding RHS components. Magnetic histories are not stored.

The packed tensor contains a contiguous signed-INT8 payload followed by a four-byte-aligned contiguous FP32 scale array. Backward maps one CUDA block to one compression tile, loads each E/R scale once into shared memory, decodes each value in a register, and immediately accumulates the epsilon/conductivity gradient. It does not allocate or write a reconstructed global-memory history. use_async_offload=True, CPU execution, and a non-float32 wavefield_storage_dtype are explicitly incompatible with "int8".

The current shared-memory maximum reduction remains available as int8_reduction_backend="current". The RTX 4090 audit selected CUB BlockReduce for the default 64-voxel tiles; "auto" falls back to the current tree for other supported power-of-two tile volumes.

E_saved is an opaque one-dimensional torch.int8 packed tensor in this mode. For diagnostics only, reconstruct it with DeepGPR.decompress_wavefield_history(E_saved, original_shape, block_size). This helper materializes FP32 and is never called by autograd backward.

CUDA Backend Build

Build the CUDA shared library from the repository root on Linux. -lineinfo preserves source correlation for Nsight; -Xptxas=-v prints registers and spill stores/loads for each kernel.

nvcc -std=c++14 -O3 -lineinfo -Xptxas=-v -arch=sm_89 --shared -Xcompiler -fPIC \
  -o src/DeepGPR/lib/deepgpr.so src/DeepGPR/lib/deepgpr.cu

The command above is the audited RTX 4090 build path. Select the matching architecture when building for a different GPU; do not add --use_fast_math unless a separate numerical and gradient validation explicitly permits it.

The new INT8 path deliberately keeps native ABI 6 because the forward/backward C signatures are unchanged. A capability symbol prevents an older ABI-6 CUDA library from accepting the packed storage code and corrupting memory. External PML also requires the deepgpr_supports_external_pml capability on CPU and CUDA, including the INT8 gradient path. Rebuild native libraries on each target platform; older binaries are rejected for PML runs rather than silently dropping physical edge gradients.

CPU Backend Build

The CPU backend is a plain C shared library and is built with OpenMP by default. Build it into src/DeepGPR/lib before running with device='cpu'. You can control CPU thread count with OMP_NUM_THREADS.

# Linux
cc -std=c99 -O3 -fopenmp -fPIC -shared -o src/DeepGPR/lib/deepgpr_cpu.so src/DeepGPR/lib/deepgpr_cpu.c

# macOS
brew install libomp
LIBOMP_PREFIX="$(brew --prefix libomp)"
LIBOMP_RUNTIME_NAME="/opt/llvm-openmp/lib/libomp.dylib"
cc -std=c99 -O3 -Xpreprocessor -fopenmp -DDEEPGPR_USE_OPENMP -I"$LIBOMP_PREFIX/include" -L"$LIBOMP_PREFIX/lib" -fPIC -shared -o src/DeepGPR/lib/deepgpr_cpu.dylib src/DeepGPR/lib/deepgpr_cpu.c -lomp
cp "$LIBOMP_PREFIX/lib/libomp.dylib" src/DeepGPR/lib/libomp.dylib
install_name_tool -id "$LIBOMP_RUNTIME_NAME" src/DeepGPR/lib/libomp.dylib
LIBOMP_DEP="$(otool -L src/DeepGPR/lib/deepgpr_cpu.dylib | awk '/libomp\.dylib/ {print $1; exit}')"
install_name_tool -change "$LIBOMP_DEP" "$LIBOMP_RUNTIME_NAME" src/DeepGPR/lib/deepgpr_cpu.dylib
codesign --force --sign - src/DeepGPR/lib/libomp.dylib
codesign --force --sign - src/DeepGPR/lib/deepgpr_cpu.dylib

On Windows, build src\DeepGPR\lib\deepgpr_cpu.dll with MSVC:

cl /LD /O2 /openmp /Fe:src\DeepGPR\lib\deepgpr_cpu.dll src\DeepGPR\lib\deepgpr_cpu.c

5. Field Variable States (Checkpoints / Initial Fields)

For starting a forward simulation from scratch ($t=0$), these three parameters should be passed as None (the system will automatically initialize zero-tensors).

checkpoint_initial_field accepts the same unextended physical model and PML settings and allocates matching full-grid states. Pass returned E, H, and all 24 PML tensors back unchanged in shape. The solver only extends materials; it never pads checkpoint states again. Continuation requires identical model geometry, PML settings, and shot ordering, with consistent materials, grid spacing, and time step for that trajectory.

The native solver advances input states in place. When a state is retained for PyTorch checkpoint recomputation or reused in another branch, clone every tensor before passing it to compute:

# Inside a checkpointed segment; boundary contains E, H, then all 24 PML states.
work = tuple(t.clone() for t in boundary)
result = DeepGPR.compute(
    device=device, dx=dx, dt=dt,
    eps_r=eps_r, sigma=sigma,  # physical model on every segment
    source_amplitudes=segment_source,
    source_location=source_location, receiver_location=receiver_location,
    pmlthick=pmlthick, E=work[:3], H=work[3:6], PML=work[6:],
)
return (*result[1], *result[2], *result[3], result[-1])

Keep material values fixed throughout one segmented trajectory and its backward pass. After the FWI optimizer step, start the next forward simulation from zero or from an initial state appropriate to that new model. See the checkpoint example.

📤 Return Values

The function returns a tuple of 5 elements. These are used to extract synthetic data, initiate the gradient flow for backpropagation, or serve as initial parameters (E, H, PML) for subsequent time-stepped calculations.

return E_saved, (Ex, Ey, Ez), (Hx, Hy, Hz), (x0EPhi1...zmHPhi2), receiver_amplitudes

1. E_saved: The pre-update electric field history E^n saved for gradient calculation and diagnostics. An internal R_saved tensor stores the corresponding discrete right-hand side R^n.
   • With save_wavefield_history=False, E_saved is a zero-length tensor and neither E nor R history storage/compression kernels are launched.
   • Shape when mode=2: (nt_saved, nstep, Nx, Ny, Nz), storing Ez only.
   • Shape when mode=3: (3, nt_saved, nstep, Nx, Ny, Nz), storing components in [Ex, Ey, Ez] order.
   • Nx, Ny, Nz include external PML. Histories, saved files, and full E/H/PML states keep that grid; only material gradients are cropped. To plot a physical history, slice spatial axes with [px0:px0+nx, py0:py0+ny, pz0:pz0+nz]. Memory estimates use the extended grid.
   • nt_saved depends on nt and model_gradient_sampling_interval.
   • Dtype is selected by wavefield_storage_dtype when compression is disabled. With wavefield_compression="int8", this is an opaque packed one-dimensional torch.int8 tensor containing values and FP32 scales.
   • Set save_forward_wavefield_path="/path/to/output" to save a CPU-loadable .pt file. Uncompressed modes save the tensor directly. INT8 mode saves a dictionary containing wavefield, compression, block_size, and uncompressed_shape so diagnostics can reconstruct it safely.
2. (Ex, Ey, Ez): The 3D electric field state at the final time step.
3. (Hx, Hy, Hz): The 3D magnetic field state at the final time step.
4. (PML_Tuple): A tuple of 24 Tensors recording the final time step state of the PML auxiliary $\Phi$ variables.
5. receiver_amplitudes: The core output. The waveform signals recorded by the receivers over the entire simulation time.
   • Shape: (nstep, nt, nrx)
   • Meaning: [Shot Index, Time Step, Receiver Index]. This output is sliced to the component specified by receiver_component (or its deprecated alias reciever_direction).

Tests and Verification

Fast unit tests and the complete numerical verification notebooks are kept in the single tests directory:

python -m unittest discover -s tests -p "test_*.py"
python tests/run_notebook.py tests/00_local_backend_and_contracts.ipynb

Run notebooks 00 through 09 in numeric order, followed by 99_verification_summary.ipynb. See tests/README.md for the complete matrix and CUDA verification options.

Cite information

If you find our codes useful, please kindly cite this article. Thanks.

@article{liu2026fast,

title={Fast ground penetrating radar dual-parameter full waveform inversion method accelerated by hybrid compilation of CUDA kernel function and PyTorch},

author={Liu, Lei and Song, Chao and He, Liangsheng and Wang, Silin and Feng, Xuan and Liu, Cai}, journal={Computers & Geosciences},

pages={106101},

year={2026},

publisher={Elsevier}

}