SPECFEM3D: 3D Elastic Wave Simulation

• SPECFEM3D is an open-source high-order spectral element method framework for 3D elastic wave propagation, using matrix-free discretization and GLL nodes.
• It simulates complex geological structures with tensor-product operations on hexahedral meshes and explicit time marching, ensuring accurate free-surface and interface handling.
• Performance is boosted through SME-aware batched matrix kernels and graph coloring, achieving significant speedups on modern HPC architectures.

Searching arXiv for recent and foundational SPECFEM3D-related papers to ground the article. SPECFEM3D is an open-source, matrix-free, high-order spectral element method (SEM) framework for 3D elastic wave propagation on hexahedral meshes, widely used in seismology and HPC. Within the SPECFEM family originally developed by Komatitsch and Tromp, it is associated with weak-form elastodynamics, tensor-product Gauss–Lobatto–Legendre (GLL) discretization, diagonal mass matrices, explicit time marching, and mesh-based representation of complex interfaces and topography (Wang et al., 11 Jun 2026, Gharti et al., 2017).

1. Numerical formulation

In the SEM formulation associated with SPECFEM3D, the governing physics is the elastodynamic wave equation for the displacement field $\mathbf{u}$ in heterogeneous media,

$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$

with the isotropic constitutive relation

$$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$

The weak form is discretized by high-order tensor-product Lagrange polynomials at GLL nodes, and the nodal GLL quadrature yields a diagonal mass matrix, which makes explicit time marching efficient. In this setting, free-surface boundary conditions are naturally satisfied, and unstructured meshes allow complex interfaces and topography to be represented accurately (Gharti et al., 2017).

A complementary description emphasizes the tensor-product structure on hexahedral elements. For polynomial degree $N$, the nodal basis is

$$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$

with interpolation nodes

$$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$

Because of the Kronecker-delta property $l_i(\xi_j)=\delta_{ij}$, directional derivatives factor into 1D operations. At the nodal points this reduces the element-local work to sequences of small dense matrix operations. For the elastic wave equation,

$$\rho\,\partial_t^2 \mathbf{s} = \nabla\cdot \mathbf{T} + \mathbf{f},\qquad \mathbf{T}=\mathbf{C}:\nabla \mathbf{s},$$

the stiffness-operator evaluation requires 3 batched small matrix multiplications to compute spatial derivatives and stresses, then another 3 batched small matrix multiplications to apply the weak form. The mass matrix is diagonal to machine precision due to GLL quadrature, so the major cost is the matrix-free stiffness operator (Wang et al., 11 Jun 2026).

2. Meshes, geometry, and boundary treatment

A defining practical feature of SPECFEM3D is its use of hexahedral meshes for realistic geological and engineering geometries. The method is presented as especially well suited to 2D and 3D heterogeneous media because material parameters are assigned elementwise and can vary sharply across the mesh. In the examples discussed in the literature, this capability is used for underground mines with air, host rock, ore, and voids; for real topography on unstable slopes; and for laboratory-scale cylindrical specimens with weak anisotropy (Gharti et al., 2017).

Topography is a central use case. For the Åknes rock slope, the free surface follows a Digital Elevation Map, and the SEM automatically satisfies the free-surface condition on this irregular boundary. The same source stresses that finite-difference methods require specialized boundary approximations and suffer accuracy loss for rough topography. In this sense, SPECFEM3D’s geometric strength lies not only in resolving internal heterogeneity but also in preserving boundary fidelity at complex free surfaces (Gharti et al., 2017).

The principal workflow limitation identified for SPECFEM3D is hexahedral meshing. In the Pyhäsalmi mine example, full 3D hexahedral meshing of the complete mine is described as impractical; the model is therefore simplified to include only two major stopes, and the geometry is decomposed into 78 volumes before meshing in CUBIT. The resulting mesh has 107,712 spectral elements and 7,161,572 spectral nodes, and it is partitioned into 24 domains for parallel processing using SCOTCH. This is presented as a typical SPECFEM3D workflow issue: automatic hexahedral meshing of complicated mine geometries is difficult, and small cavities and inclusions are usually not worth honoring if they would explode the mesh size (Gharti et al., 2017).

A related practical result concerns voids and air-filled cavities. In the 2D Pyhäsalmi model, including air in the stopes forces a very fine mesh because air has $V_s = 0$ and very low $V_p$. Two meshes were tested: one including air in the stopes, with average element size 2 m and 53,071 spectral elements, and one excluding the stopes and air cavities, with 40,594 spectral elements. Both simulations were stable, and the synthetic wavefields were very similar, with only negligible waveform differences. The conclusion drawn there is that, if acoustic propagation inside the voids is not the target, removing the cavities is a safe and cost-effective choice, reducing the element count by about 30% in that 2D example (Gharti et al., 2017).

3. Coupled domains and interface treatments

Although SPECFEM3D is primarily described here as a 3D elastic SEM code, its scientific neighborhood includes fluid-solid and elastic-acoustic coupling. A closely related formulation is the 3D coupled elastic–acoustic wave problem in a symmetric second-order form, with displacement $$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$0 in the elastic domain and velocity potential $$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$1 in the acoustic domain, defined by

$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$2

The strong form uses

$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$3

$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$4

with interface conditions

$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$5

In the DG Spectral Element variant, the interface coupling is handled so that elastic and acoustic meshes may be non-matching at $$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$6, whereas classic conforming SEM implementations such as the original SPECFEM-style fluid-solid coupling often assume matching interface grids or use specialized interface treatments (Antonietti et al., 2019).

This comparison clarifies a recurrent distinction. SPECFEM3D is repeatedly grouped with high-order 3D wave-propagation solvers for realistic Earth models, but it is not numerically interchangeable with all of them. In the fully coupled Palu earthquake–tsunami study, SPECFEM3D is used as a conceptual reference point for high-order 3D wave propagation in realistic geometry, whereas the actual solver is SeisSol, implemented with ADER-DG on adaptive unstructured tetrahedral meshes. There, SPECFEM3D is associated with spectral elements on hexahedral meshes and global or structured-element-style assembly workflows, while SeisSol uses face-based Riemann solvers, multi-rate local time stepping, and a coupled treatment of dynamic rupture, elastic wave propagation in the solid Earth, acoustic wave propagation in the ocean, and tsunami generation with gravity (Krenz et al., 2021).

The consequence is methodological rather than hierarchical. A plausible implication is that SPECFEM3D occupies a central position in the spectral-element branch of large-scale wave simulation, while adjacent frameworks trade that structure for different interface flexibility, coupling mechanisms, or mesh types.

4. Applications across scales

SPECFEM3D has been used to model full seismic and acoustic wavefields in microearthquake environments with strong small-scale structural and material heterogeneities. The documented examples span an underground ore mine, a steep topographic slope, and a laboratory specimen under triaxial loading. These cases are not introduced as algorithmic novelties; rather, they demonstrate the range of environments in which the standard SPECFEM workflow remains viable (Gharti et al., 2017).

In the Pyhäsalmi underground ore mine in Finland, the volcanogenic massive sulphide deposit extends to about 1.4 km depth, and the in-mine seismic network consists of 18 geophones, including 6 three-component stations, with sampling rate as high as 3 kHz. The 3D simulation of a real observed microseismic event uses a simplified geometry with two major stopes excluded from meshing. The source is represented by a Ricker wavelet at 200 Hz and a full moment tensor with both isotropic and deviatoric parts. Comparison with observed waveforms shows generally good agreement in first-arrival times and coda durations, but not a precise match, which is attributed to the complex medium and uncertain source mechanism. The wavefield exhibits shadow zones caused by voids and source radiation pattern, and the analysis stresses the importance of 3-component recordings for reliable interpretation (Gharti et al., 2017).

At the Åknes unstable rock slope in western Norway, the unstable mass is about 35–40 million m$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$7, with motion of roughly 6 cm/year on average and up to 14 cm/year, in a setting that poses a tsunami hazard. Because no reliable velocity model is available, the model isolates the effect of true topography by using a homogeneous medium with

$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$8

The source is an explosion moment tensor on the free surface with a Ricker wavelet of central frequency 120 Hz. The mesh has average element size 10 m, 107,712 spectral elements, 7,161,572 spectral nodes, and is partitioned into 8 domains. The resulting wavefield is strongly complicated by multiple reflections and conversions at the topography; S waves are generated by free-surface interactions and can be stronger than the P waves (Gharti et al., 2017).

At laboratory scale, the same machinery is applied to a Vosges sandstone cylinder of diameter 25.4 mm and height 63.5 mm, with a central borehole of diameter 5.2 mm. The triaxial experiment produces about 2500 detected acoustic-emission events, recorded by 12 piezoelectric sensors at 10 MHz sampling rate. The homogeneous model uses

$$\rho \,\ddot{\mathbf{u}} = \nabla \cdot \boldsymbol{\sigma} + \mathbf{f},$$9

and weak anisotropy is introduced with Thomsen parameters

$$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$0

The mesh contains 82,240 spectral elements and 5,462,768 spectral nodes, with average element size 1.5 mm. Synthetic isotropic and anisotropic waveforms broadly match observations; anisotropy produces only small delays, and shear-wave splitting is not clearly visible, indicating that the anisotropy is weak (Gharti et al., 2017).

5. Performance and architecture-specific optimization

SPECFEM3D is also treated as a representative real-world HPC workload for SEM. On the LX2 ARMv9-A processor, the target architecture is described as having 304 cores total, 4 NUMA domains per die, 4 GB HBM per NUMA node, shared DDR, and SME support through ZA tiles. The architectural challenge is that SEM tensor-product kernels consist of small matrix multiplications; for practical polynomial orders, the matrices are often too small to saturate the SME tile efficiently, so a naïve SME port would underutilize the hardware (Wang et al., 11 Jun 2026).

The proposed optimization strategy has three parts. First, an SME-aware batched small-matrix kernel redesigns the tensor-product operators so that the $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$1-, $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$2-, and $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$3-direction derivatives better match SME’s outer-product execution model. Second, a memory-aware hybrid MPI+OpenMP execution model maps MPI processes to NUMA domains and uses OpenMP threads within a NUMA domain in order to reduce duplicated data structures per process and avoid exhausting the 4 GB HBM per NUMA node. Third, a dispersion-based iso-accuracy study of the $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$4 tradeoff is used to reassess the practical discretization choice rather than only the kernel micro-optimization (Wang et al., 11 Jun 2026).

The same study reports that SPECFEM3D traditionally uses atomic updates for shared nodal values, but atomics are expensive on CPUs because they rely on LL/SC or CAS retry loops and disrupt pipelining. The replacement is graph coloring: build an element adjacency graph from mesh connectivity, color elements so that elements in the same color class do not share nodes, then execute color by color so that no atomics are needed within each color. Inter-process communication is overlapped by a dedicated communication/progress thread handling nonblocking MPI progress, with one CPU core per MPI process reserved for progress. The benefit is described as real but limited in practice by load imbalance and very small halo sizes (Wang et al., 11 Jun 2026).

The reported end-to-end performance gains are strongly polynomial-order dependent. Compared with the original SPECFEM3D, vectorization alone improves performance by $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$5 at $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$6, $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$7 at $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$8, and $$\boldsymbol{\sigma} = \lambda (\nabla \cdot \mathbf{u}) \mathbf{I} + 2\mu \boldsymbol{\varepsilon}(\mathbf{u}), \qquad \boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}\left(\nabla \mathbf{u} + \nabla \mathbf{u}^{T}\right).$$9 at $N$0. Graph coloring adds $N$1 at $N$2, $N$3 at $N$4, and has a slightly negative impact at $N$5 ($N$6). Adding SME on top yields stiffness-operator speedups of $N$7 at $N$8, $N$9 at $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$0, and $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$1 at $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$2, translating into full-application speedups of $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$3 at $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$4, $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$5 at $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$6, and $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$7 at $$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$8. The abstract summarizes the overall full-application gain as 4–6$$l_{ijk}(\xi,\eta,\zeta)=l_i(\xi)\,l_j(\eta)\,l_k(\zeta),$$9 at fixed polynomial order (Wang et al., 11 Jun 2026).

The iso-accuracy study is equally consequential. With waveform error targeted below 1% after 8 s and a fixed time step

$$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$0

the maximum admissible element sizes are reported as 20 m for $$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$1, 50 m for $$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$2, and 130 m for $$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$3. Moving from $$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$4 to $$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$5 gives 3.18$$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$6 less memory and 3.85$$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$7 shorter time-to-solution; moving from $$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$8 to $$\boldsymbol{x}_{ijk}=(\xi_i,\xi_j,\xi_k),\qquad i,j,k=0,\dots,N.$$9 gives another 1.8$l_i(\xi_j)=\delta_{ij}$0 less memory and another 1.4$l_i(\xi_j)=\delta_{ij}$1 shorter time-to-solution. The stated interpretation is that SME changes not only kernel efficiency but also the performance-favorable operating point along the SEM $l_i(\xi_j)=\delta_{ij}$2 iso-accuracy frontier, shifting it toward higher polynomial orders (Wang et al., 11 Jun 2026).

6. Position within high-order wave simulation

SPECFEM3D belongs to a broader ecosystem of high-order 3D wave-propagation solvers, but the literature repeatedly distinguishes it from neighboring methods rather than treating them as interchangeable. In relation to DG spectral element methods, the important common ground is high-order tensor-product polynomial bases on hexahedra, explicit time integration, 3D wave propagation, and absorbing boundary conditions. The significant difference is interface treatment: the DGSE approach localizes discontinuities at selected interfaces and allows independent meshing of elastic and acoustic subdomains, whereas SPECFEM3D is presented as a conforming spectral-element framework with more constrained interface handling (Antonietti et al., 2019).

In relation to SeisSol, the distinction is sharper. The fully coupled Palu study explicitly notes that it is not a SPECFEM3D paper, even though it addresses the same broad class of problems—3D seismic wave propagation in realistic geometry. There, SPECFEM3D is the reference point for spectral elements, while SeisSol is the vehicle for a more difficult multiphysics problem: a 3D fully coupled earthquake–tsunami system involving non-linear dynamic rupture, elastic wave propagation in the Earth, acoustic wave propagation in the ocean, and gravity-modified free-surface dynamics. The comparison is therefore not a replacement claim; it is a statement that extensions to compressible ocean acoustics, tsunami gravity waves, and dynamic rupture may motivate a different discretization and solver infrastructure (Krenz et al., 2021).

Several common misconceptions are corrected by these comparisons. One is that all high-order wave solvers for realistic geology are numerically equivalent; the cited work instead separates spectral elements on hexahedral meshes from ADER-DG on tetrahedra and from DGSE interface formulations. Another is that lower polynomial order is automatically the safe performance choice; on ARM multicore CPUs with SME, the reported iso-accuracy results favor higher $l_i(\xi_j)=\delta_{ij}$3 because they reduce both working-set size and time-to-solution at fixed accuracy. A third is that every geometrical feature must be explicitly meshed; in the mine example, excluding stopes and air cavities is shown to preserve the external wavefield closely while reducing cost when acoustic propagation inside the voids is not of interest (Wang et al., 11 Jun 2026, Gharti et al., 2017).

Taken together, these accounts place SPECFEM3D as a canonical high-order SEM code for 3D elastic wave propagation in complex geometry, with a mature workflow centered on hexahedral meshing, GLL-based diagonal-mass discretization, explicit time marching, and large-scale parallel execution. Its limitations and adjacent alternatives are not peripheral details; they define the methodological boundaries within which SPECFEM3D is most effective.