I’m currently working on a GPU-accelerated 2D elastic wave simulation code in Julia,
based on staggered-grid finite-difference discretization.
The original motivation is seismic forward modeling, but I’m also interested in
making the implementation reasonably idiomatic and reusable in the Julia ecosystem.
At the moment, the project supports:
2D elastic wave equation on a staggered grid
Hybrid absorbing boundary conditions for elastic waves
Explicit time stepping
Explicit CUDA.jl-based GPU kernels for field updates
The code is currently research-oriented:
it runs end-to-end on the GPU and produces reasonable results,
but the public API and package structure are still evolving.
The repository is here:
This is very much a work in progress.
Possible future directions
Some possible future directions I’m considering (in no particular commitment order):
Support for more general free-surface geometries
(e.g. irregular topography beyond simple flat boundaries)
Extension from 2D to 3D elastic wave modeling,
once the core abstractions are stable enough
Multi-GPU execution via domain decomposition,
potentially using MPI and CUDA-aware communication
I’d really appreciate feedback on a few specific points:
API design
Does the current separation between configuration structs and simulation routines
look idiomatic from a Julia perspective?
GPU kernel organization
Is it reasonable to keep explicit CUDA kernels for the core updates,
or would you recommend a higher-level abstraction (e.g. KernelAbstractions.jl)
at this stage?
Project scope
From a community perspective, does this look more like:
a reusable library, or
an application / research code that should stay more opinionated?
Related work
Are there existing Julia packages or projects I should study or align with?
Thanks for your time, and any suggestions are very welcome.
I wonder, however, whether for-loops as in the function below would benefit from the use of Tullio.jl?
function update_v_core!(vx, vz, txx, tzz, txz, rho_vx, rho_vz, a, nx, nz, dtx, dtz, M_order)
@tturbo for j in (M_order+1):(nz-M_order)
for i in (M_order+1):(nx-M_order)
dtxxdx, dtxzdz, dtxzdx, dtzzdz = 0.0f0, 0.0f0, 0.0f0, 0.0f0
# High-order staggered finite difference stencil
for l in 1:M_order
# Update vx (located at i, j)
# Derivatives of txx (i+0.5, j) and txz (i, j+0.5)
dtxxdx += a[l] * (txx[i+l-1, j] - txx[i-l, j])
dtxzdz += a[l] * (txz[i, j+l-1] - txz[i, j-l])
...
end
...
end
end
end
Thank you for the suggestion!
I haven’t tried Tullio.jl yet, but it sounds promising — the Einstein notation would make the stencil code much cleaner.I will try it.
I have other naughty and possibly irrelevant questions.
Q1: You currently hand code the explicit time integration with a fixed time step. Is it relevant here to exploit a variable time step? Can functions already available be turned into a right-hand function for DifferentialEquations.jl ? What is required to define a minimal working example (MWE) in this direction?
Q2: Does your application benefit from implicit time-integration (thus requiring a linear solve at each time step)?
Q3: Is your application limited to linear wave propagation or do you instead have tackling non-linearities (non-elastic) in mind (thus requiring a non-linear solve at each time step)?
Q4: have you considered the possibility of coupling the physics of wave propagation to other physical phenomena, e.g., the transport of a passive scalar with the wave field? [In our application we have shalllow water flow coupled with sediment transport].
Thank you for the thoughtful questions and for taking the time to review the project.
Fomo is mainly designed to generate seismic data for large-scale models, with a strong focus on performance. Because of CFL constraints and the nature of wave propagation, I use explicit time stepping with a fixed time step. This is a common approach in geophysics, and variable time stepping usually does not bring much benefit here. Using DifferentialEquations.jl is possible in principle, but it is not a goal of the project at the moment.
In the future, I plan to make the code more modular, so that different wave equations and options can be selected more easily (e.g. acoustic vs elastic, isotropic vs anisotropic media, different boundary conditions such as free surfaces).
For now, the scope is limited to linear wave propagation. Implicit methods, nonlinear effects, and multi-physics coupling are not considered yet. The current goal is to make Fomo a simple, efficient, and reliable seismic forward-modeling tool.
I’ve now turned the project into a Julia package. Developing this in Julia has been very convenient, especially for structuring the code and experimenting with different ideas.
Are the 2D elastic models only defined as property grids?
Can we include any type of layer interfaces, as they control how waves are reflected, transmitted, and refracted?
These interfaces are important in some cases for the accuracy of the simulation results.
Thanks for sharing this package.
I haven’t considered a velocity model generation approach based on explicit interfaces for now; in fact, this is essentially a 2D interpolation problem.
If you want to perform forward modeling based on interfaces, you could use SPECFEM2D or SPECFEM3D. These are spectral-element–based simulation tools and are very well suited for interface-based forward modeling.
Thank you for taking the time to look through my code.