Using a Gmsh mesh in SPECFEM++ (3D)¶
As with any finite element solver, SPECFEM++ requires a mesh generation step to
discretize the domain. The spectral element solver requires a hexahedral mesh and
does not support tetrahedral meshes. In this example we use
Gmsh — a free, cross-platform mesh generator — to create a
simple 3D box domain, export it using a Python script, and run the simulation using
the xdecompose_mesh partitioner followed by specfem3d.
Note
This cookbook is not an in-depth Gmsh tutorial. Its goal is to show how to define a Gmsh mesh so that SPECFEM++ can read it. For more information on meshing with Gmsh please refer to the Gmsh tutorials.
Note
All cookbook files can be copy and pasted from the code blocks below, or you can download a zip file containing all the files needed to run this example.
Download GMSH 3D cookbook here
Or download using command line:
# Using curl:
curl -O https://specfem2d-kokkos.readthedocs.io/en/devel/sections/cookbooks/wavepropagation/dim3/Gmsh/gmsh_3d_cookbook.zip
# Using wget:
wget https://specfem2d-kokkos.readthedocs.io/en/devel/sections/cookbooks/wavepropagation/dim3/Gmsh/gmsh_3d_cookbook.zip
Setting up your workspace¶
Install Gmsh and its Python API (if not already available):
pip install gmsh
Create a working directory and make sure the SPECFEM++ executables are on your
PATH.
mkdir -p ~/specfempp-examples/gmsh-halfspace-3d
cd ~/specfempp-examples/gmsh-halfspace-3d
which specfem3d
If the command does not return a path, add the executable directory to PATH:
export PATH=$PATH:<PATH TO SPECFEM++ DIRECTORY>/bin
Create the output directories and placeholder files:
mkdir -p MESH
mkdir -p DATABASES_MPI
mkdir -p OUTPUT_FILES/results
mkdir -p DATA
touch specfem_config.yaml
touch sources.yaml
Creating the geometry¶
We will create a 10 km × 10 km × 5 km box domain with the free surface at z = 0 and the bottom of the model at z = -5000 m.
Save the following as create_mesh.py in your working directory:
import gmsh
gmsh.initialize()
gmsh.model.add("halfspace_3d")
## Create a 10 km x 10 km x 5 km box: x:[0,10000], y:[0,10000], z:[-5000,0]
box = gmsh.model.occ.addBox(0, 0, -5000, 10000, 10000, 5000)
gmsh.model.occ.synchronize()
## ---------------------------------------------------------------
## Physical groups
## Volume group for the elastic material
## ---------------------------------------------------------------
gmsh.model.addPhysicalGroup(3, [box], name="material_1")
## ---------------------------------------------------------------
## Surface groups for boundary conditions
## Use the bounding box of each surface to identify the face.
## ---------------------------------------------------------------
eps = 1.0
for dim, tag in gmsh.model.getEntities(2):
xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(dim, tag)
if xmax < eps:
gmsh.model.addPhysicalGroup(2, [tag], name="abs_xmin")
elif xmin > 10000 - eps:
gmsh.model.addPhysicalGroup(2, [tag], name="abs_xmax")
elif ymax < eps:
gmsh.model.addPhysicalGroup(2, [tag], name="abs_ymin")
elif ymin > 10000 - eps:
gmsh.model.addPhysicalGroup(2, [tag], name="abs_ymax")
elif zmax < -5000 + eps:
gmsh.model.addPhysicalGroup(2, [tag], name="abs_bottom")
else: ## z = 0: free surface
gmsh.model.addPhysicalGroup(2, [tag], name="free_or_abs_zmax")
## ---------------------------------------------------------------
## Structured transfinite hex mesh at 500 m element size
## ---------------------------------------------------------------
## Number of nodes along each edge (nodes = elements + 1)
n_horiz = 21 ## 20 elements across 10 km = 500 m
n_vert = 11 ## 10 elements across 5 km = 500 m
for dim, tag in gmsh.model.getEntities(1):
xmin, ymin, zmin, xmax, ymax, zmax = gmsh.model.getBoundingBox(dim, tag)
length = max(abs(xmax - xmin), abs(ymax - ymin), abs(zmax - zmin))
n = round(length / 500) + 1
gmsh.model.mesh.setTransfiniteCurve(tag, n)
for dim, tag in gmsh.model.getEntities(2):
gmsh.model.mesh.setTransfiniteSurface(tag)
gmsh.model.mesh.setRecombine(dim, tag)
gmsh.model.mesh.setTransfiniteVolume(box)
gmsh.model.mesh.generate(3)
gmsh.write("halfspace.msh")
gmsh.finalize()
print("Mesh written to halfspace.msh")
Run the script to generate the mesh file:
python create_mesh.py
After the script completes you should see halfspace.msh in the working
directory. You can open it in the Gmsh GUI to inspect the mesh:
gmsh halfspace.msh
Note
The transfinite structured meshing algorithm requires that opposite faces of
the volume are meshed with the same topology. For a simple box this is
guaranteed automatically. For more complex geometries you may need to specify
the source and target surfaces for the transfinite sweep manually using
gmsh.model.mesh.setTransfiniteVolume(tag, corner_points).
After meshing you should see approximately 4 000 hexahedral elements (20 × 20 × 10).
Understanding physical groups¶
Physical groups in Gmsh serve the same role as block names in CUBIT: they tag mesh entities so the export script knows how to route them to the correct output files. The naming convention mirrors the CUBIT workflow exactly:
Group name (dim) |
Output file |
|---|---|
|
element connectivity, material mapping, velocity file |
|
|
|
|
|
|
|
|
|
|
|
|
Tip
The bounding-box approach in create_mesh.py automatically identifies
each face regardless of the internal surface tag numbers Gmsh assigns.
This makes the script robust to geometry modifications.
Exporting the mesh¶
The export script is available in the SPECFEM++ repository at
<PATH TO SPECFEM++ DIRECTORY>/scripts/export_gmsh3d.py.
Unlike the CUBIT export script, export_gmsh3d.py is a standard Python script
that reads the saved .msh file — it does not need to run inside Gmsh.
Add it to your path or copy it to your working directory, then run:
import sys
sys.path.insert(0, "<PATH TO SPECFEM++ DIRECTORY>/scripts")
from export_gmsh3d import export2SPECFEM3D_gmsh, MaterialSpec
## Material properties for material_1 (elastic granite)
materials = {
1: MaterialSpec(
mat_id = 1,
domain_id = 2, ## 2 = elastic
rho = 2700.0, ## density (kg/m^3)
vp = 6000.0, ## P-wave velocity (m/s)
vs = 3500.0, ## S-wave velocity (m/s)
q_kappa = 9999.0, ## no attenuation
q_mu = 9999.0,
),
}
export2SPECFEM3D_gmsh("halfspace.msh", materials, outdir="MESH")
python export.py
Alternatively, use the command-line interface directly:
python <PATH TO SPECFEM++ DIRECTORY>/scripts/export_gmsh3d.py \
halfspace.msh MESH \
--mat 1 2 2700 6000 3500
This will print status messages and generate the following files inside the MESH
directory:
nodes_coords_file— node coordinatesmesh_file— hexahedral element connectivitymaterials_file— element-to-material mappingnummaterial_velocity_file— material velocity propertiesabsorbing_surface_file_xmin— xmin absorbing face elementsabsorbing_surface_file_xmax— xmax absorbing face elementsabsorbing_surface_file_ymin— ymin absorbing face elementsabsorbing_surface_file_ymax— ymax absorbing face elementsabsorbing_surface_file_bottom— bottom absorbing face elementsfree_or_absorbing_surface_file_zmax— top free surface elements
Note
The domain_id in MaterialSpec controls the physics:
1 = acoustic (fluid), 2 = elastic (solid).
Set vs = 0 for acoustic materials.
Running xdecompose_mesh¶
xdecompose_mesh partitions the mesh text files into a binary database that
specfem3d can read. We first define the decomposer parameters in a Par_file.
#-----------------------------------------------------------
#
# Parameters for xdecompose_mesh
#
#-----------------------------------------------------------
# number of MPI partitions to create
NPROC = 1
# number of nodes per hexahedral element (8 or 27)
NGNOD = 8
# partitioner: 1=SCOTCH (default), 2=METIS, 3=PATOH, 4=ROWS_PART
PARTITIONING_TYPE = 1
# local-time stepping (must match solver Par_file)
LTS_MODE = .false.
# attenuation (affects database layout)
ATTENUATION = .false.
# C-PML absorbing boundary conditions
PML_CONDITIONS = .false.
# save VTK mesh files for quality checking
SAVE_MESH_FILES = .true.
# ADIOS database output (not supported in serial decomposer -- must be .false.)
ADIOS_FOR_DATABASES = .false.
# HDF5 database output
HDF5_ENABLED = .false.
# injection coupling
COUPLE_WITH_INJECTION_TECHNIQUE = .false.
# wavefield discontinuity interface
IS_WAVEFIELD_DISCONTINUITY = .false.
# output directory for the per-partition database files
LOCAL_PATH = ./DATABASES_MPI
#-----------------------------------------------------------
#
# Input mesh file paths
# (relative to the directory from which xdecompose_mesh is invoked)
#
#-----------------------------------------------------------
# required
NODES_COORDS_FILE = MESH/nodes_coords_file
MESH_FILE = MESH/mesh_file
MATERIALS_FILE = MESH/materials_file
NUMMATERIAL_VELOCITY_FILE = MESH/nummaterial_velocity_file
# absorbing boundary surfaces (comment out any that are not present)
ABSORBING_SURFACE_FILE_XMIN = MESH/absorbing_surface_file_xmin
ABSORBING_SURFACE_FILE_XMAX = MESH/absorbing_surface_file_xmax
ABSORBING_SURFACE_FILE_YMIN = MESH/absorbing_surface_file_ymin
ABSORBING_SURFACE_FILE_YMAX = MESH/absorbing_surface_file_ymax
ABSORBING_SURFACE_FILE_BOTTOM = MESH/absorbing_surface_file_bottom
FREE_OR_ABSORBING_SURFACE_FILE_ZMAX = MESH/free_or_absorbing_surface_file_zmax
Key parameters to note:
NPROC— number of MPI partitions. Set to 1 here for a serial run. Increase it to match the number of MPI processes you plan to use withspecfem3d.NGNOD = 8— matches our HEX8 mesh.LOCAL_PATH— directory where the binary partition databases are written.The
NODES_COORDS_FILEthroughFREE_OR_ABSORBING_SURFACE_FILE_ZMAXentries point to the files generated byexport2SPECFEM3D_gmshin the previous step.
Run the decomposer from your working directory:
xdecompose_mesh -p Par_file
After a successful run you should see a binary database file inside DATABASES_MPI:
ls DATABASES_MPI/
For NPROC = 1 this produces DATABASES_MPI/Database.bin. For NPROC > 1
a subdirectory DATABASES_MPI/Database/ is created containing one
proc_<rank>.bin file per MPI rank.
Defining receivers (stations)¶
Receivers are defined in the DATA/STATIONS file. Each line specifies one station
with the format:
STATION_NAME NETWORK_CODE X Y ELEVATION BURIAL
ST001 DB 2000.0 5000.0 0.0 0.0
ST002 DB 3000.0 5000.0 0.0 0.0
ST003 DB 4000.0 5000.0 0.0 0.0
ST004 DB 6000.0 5000.0 0.0 0.0
ST005 DB 7000.0 5000.0 0.0 0.0
ST006 DB 8000.0 5000.0 0.0 0.0
These six stations are distributed along the free surface (z = 0), displaced 1–3 km from the source in the x direction.
Defining sources¶
The source is defined in sources.yaml. We place a vertical force at the centre of
the domain, 500 m below the free surface, with a 1 Hz Ricker wavelet.
number-of-sources: 1
sources:
- force:
x : 5000.0
y : 5000.0
z : -500.0
source_surf: false
fx : 0.0
fy : 0.0
fz : 1e10
Ricker:
factor: 1.0
tshift: 0.0
f0: 1.0
Configuring the solver¶
Define the simulation parameters in specfem_config.yaml:
parameters:
header:
title: 3D Gmsh external mesh -- homogeneous elastic halfspace
description: |
Material systems : Elastic domain (1)
Interfaces : None
Sources : Force source (1), Ricker wavelet, f0=1 Hz
Boundary conditions : Absorbing (xmin, xmax, ymin, ymax, bottom), Free surface (zmax)
simulation-setup:
quadrature:
quadrature-type: GLL4
solver:
time-marching:
time-scheme:
type: Newmark
dt: 0.004
nstep: 1250
simulation-mode:
forward:
writer:
seismogram:
format: "ascii"
directory: "OUTPUT_FILES/results"
# display:
# format: VTKHDF
# directory: OUTPUT_FILES/
# field: displacement
# component: z
# simulation-field: forward
# time-interval: 5
receivers:
stations: "DATA/STATIONS"
angle: 0.0
seismogram-type:
- displacement
nstep_between_samples: 1
run-setup:
number-of-processors: 1
number-of-runs: 1
databases:
mesh-database: "DATABASES_MPI/Database.bin"
sources: "sources.yaml"
A few parameters to highlight:
dt = 0.004s andnstep = 1250gives 5 s of simulation, enough to record both the direct P and S arrivals at all stations. The CFL stability limit for this mesh (500 m elements, GLL4, vp = 6000 m/s) is approximately 0.3 × 86 m / 6000 m/s ≈ 0.0043 s, sodt = 0.004s sits comfortably within the stable range.mesh-databasepoints toDATABASES_MPI/Database.bin, which is wherexdecompose_meshwrites the single-process partition database. If you usedNPROC > 1you do not need to change this path —specfem3dautomatically appends the rank suffix.
Running the solver¶
specfem3d -p specfem_config.yaml
The solver writes ASCII seismogram files for each station and component to
OUTPUT_FILES/results/.
ls OUTPUT_FILES/results/
You should see one file per station per displacement component, e.g.:
DB.ST001.BXX.semd, DB.ST001.BXY.semd, DB.ST001.BXZ.semd.
[Optional] Visualizing the wavefield¶
To visualize the wavefield propagation SPECFEM++ must be built with VTK and HDF5
support. Uncomment the display section in specfem_config.yaml:
display:
format: VTKHDF
directory: OUTPUT_FILES/
field: displacement
component: z
simulation-field: forward
time-interval: 5
Rerun the solver, then open Paraview and load the generated VTK files from the
OUTPUT_FILES directory to visualize the wavefield propagation.