Elastic & viscoelastic wave propagation comparison¶
In this example (see benchmarks/src/dim2/homogeneous-medium-flat-topography-attenuation)
we simulate wave propagation through the same 2-dimensional homogeneous elastic medium as in
Homogeneous elastic media, but run the simulation twice: once without attenuation and once with
viscoelastic attenuation enabled. Comparing the two sets of seismograms reveals the amplitude
decay and velocity dispersion introduced by attenuation.
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 Attenuation cookbook here
Or download using command line:
# Using curl:
curl -O https://specfem2d-kokkos.readthedocs.io/en/devel/sections/cookbooks/wavepropagation/dim2/homogeneous-medium-flat-topography-attenuation/attenuation_cookbook.zip
# Using wget:
wget https://specfem2d-kokkos.readthedocs.io/en/devel/sections/cookbooks/wavepropagation/dim2/homogeneous-medium-flat-topography-attenuation/attenuation_cookbook.zip
Setting up your workspace¶
Let’s start by creating a workspace from where we can run this example.
mkdir -p ~/specfempp-examples/attenuation
cd ~/specfempp-examples/attenuation
We also need to check that the SPECFEM++ executable directory is added to the
PATH.
which specfem
If the above command returns a path to the specfem executable, then the
executable directory is added to the PATH. If not, you need to add the executable
directory to the PATH using the following command.
export PATH=$PATH:<PATH TO SPECFEM++ DIRECTORY/bin>
Note
Make sure to replace <PATH TO SPECFEM++ DIRECTORY/bin> with the
actual path to the SPECFEM++ directory on your system.
Now let’s create the necessary directories for both simulation variants.
mkdir -p OUTPUT_FILES/attenuation_off/results
mkdir -p OUTPUT_FILES/attenuation_on/results
mkdir -p OUTPUT_FILES/results
touch specfem_config_attenuation_off.yaml
touch specfem_config_attenuation_on.yaml
touch source.yaml
touch Par_File_attenuation_off
touch Par_File_attenuation_on
touch topography.dat
Generating the mesh¶
Both simulation variants use the same homogeneous medium and identical mesh, so the
mesh generation steps are the same except for the output directory. We run xmeshfem2D
separately for each variant so that the STATIONS file and database are written to the
correct subdirectory.
Parameter File (without attenuation)¶
Par_File_attenuation_off¶
#-----------------------------------------------------------
#
# Simulation input parameters
#
#-----------------------------------------------------------
# title of job
title = Elastic Simulation with point source
# parameters concerning partitioning
NPROC = 1 # number of processes
# Output folder to store mesh related files
OUTPUT_FILES = OUTPUT_FILES/attenuation_off
#-----------------------------------------------------------
#
# Mesh
#
#-----------------------------------------------------------
# Partitioning algorithm for decompose_mesh
PARTITIONING_TYPE = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
# number of control nodes per element (4 or 9)
NGNOD = 9
# location to store the mesh
database_filename = OUTPUT_FILES/attenuation_off/database.bin
#-----------------------------------------------------------
#
# Receivers
#
#-----------------------------------------------------------
# use an existing STATION file found in ./DATA or create a new one from the receiver positions below in this Par_file
use_existing_STATIONS = .false.
# number of receiver sets (i.e. number of receiver lines to create below)
nreceiversets = 2
# orientation
anglerec = 0.d0 # angle to rotate components at receivers
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
nrec = 3 # number of receivers
xdeb = 2200. # first receiver x in meters
zdeb = 2200. # first receiver z in meters
xfin = 2800. # last receiver x in meters (ignored if only one receiver)
zfin = 2200. # last receiver z in meters (ignored if only one receiver)
record_at_surface_same_vertical = .true. # receivers inside the medium or at the surface (z values are ignored if this is set to true, they are replaced with the topography height)
# second receiver set
nrec = 3 # number of receivers
xdeb = 2500. # first receiver x in meters
zdeb = 2500. # first receiver z in meters
xfin = 2500. # last receiver x in meters (ignored if only one receiver)
zfin = 1900. # last receiver z in meters (ignored if only one receiver)
record_at_surface_same_vertical = .false. # receivers inside the medium or at the surface (z values are ignored if this is set to true, they are replaced with the topography height)
# filename to store stations file
stations_filename = OUTPUT_FILES/attenuation_off/STATIONS
#-----------------------------------------------------------
#
# Velocity and density models
#
#-----------------------------------------------------------
# number of model materials
nbmodels = 1
# available material types (see user manual for more information)
# acoustic: model_number 1 rho Vp 0 0 0 QKappa 9999 0 0 0 0 0 0 (for QKappa use 9999 to ignore it)
# elastic: model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0 (for QKappa and Qmu use 9999 to ignore them)
# anisotropic: model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 0 QKappa Qmu
# anisotropic in AXISYM: model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 c22 QKappa Qmu
# poroelastic: model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu
# tomo: model_number -1 0 0 A 0 0 0 0 0 0 0 0 0 0
#
# note: When viscoelasticity or viscoacousticity is turned on,
# the Vp and Vs values that are read here are the UNRELAXED ones i.e. the values at infinite frequency
# unless the READ_VELOCITIES_AT_f0 parameter above is set to true, in which case they are the values at frequency f0.
#
# Please also note that Qmu is always equal to Qs, but Qkappa is in general not equal to Qp.
# To convert one to the other see doc/Qkappa_Qmu_versus_Qp_Qs_relationship_in_2D_plane_strain.pdf and
# utils/attenuation/conversion_from_Qkappa_Qmu_to_Qp_Qs_from_Dahlen_Tromp_959_960.f90.
1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
# external tomography file
TOMOGRAPHY_FILE = ./DATA/tomo_file.xyz
# use an external mesh created by an external meshing tool or use the internal mesher
read_external_mesh = .false.
#-----------------------------------------------------------
#
# PARAMETERS FOR EXTERNAL MESHING
#
#-----------------------------------------------------------
# data concerning mesh, when generated using third-party app (more info in README)
# (see also absorbing_conditions above)
mesh_file = ./DATA/mesh_file # file containing the mesh
nodes_coords_file = ./DATA/nodes_coords_file # file containing the nodes coordinates
materials_file = ./DATA/materials_file # file containing the material number for each element
free_surface_file = ./DATA/free_surface_file # file containing the free surface
axial_elements_file = ./DATA/axial_elements_file # file containing the axial elements if AXISYM is true
absorbing_surface_file = ./DATA/absorbing_surface_file # file containing the absorbing surface
acoustic_forcing_surface_file = ./DATA/MSH/Surf_acforcing_Bottom_enforcing_mesh # file containing the acoustic forcing surface
absorbing_cpml_file = ./DATA/absorbing_cpml_file # file containing the CPML element numbers
tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
#-----------------------------------------------------------
#
# PARAMETERS FOR INTERNAL MESHING
#
#-----------------------------------------------------------
# file containing interfaces for internal mesh
interfacesfile = topography.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
xmax = 4000.d0 # abscissa of right side of the model
nx = 80 # number of elements along X
STACEY_ABSORBING_CONDITIONS = .false.
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
absorbright = .true.
absorbtop = .true.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral-element mesh
nbregions = 1 # then set below the different regions and model number for each region
# format of each line: nxmin nxmax nzmin nzmax material_number
1 80 1 60 1
#-----------------------------------------------------------
#
# DISPLAY PARAMETERS
#
#-----------------------------------------------------------
# meshing output
output_grid_Gnuplot = .false. # generate a GNUPLOT file containing the grid, and a script to plot it
output_grid_ASCII = .false. # dump the grid in an ASCII text file consisting of a set of X,Y,Z points or not
The key parameter is the material definition on line 90: QKappa=9999 Qmu=9999 — setting
both quality factors to 9999 tells the solver to treat the medium as purely elastic (no
attenuation).
Parameter File (with attenuation)¶
Par_File_attenuation_on¶
#-----------------------------------------------------------
#
# Simulation input parameters
#
#-----------------------------------------------------------
# title of job
title = Elastic Simulation with point source
# parameters concerning partitioning
NPROC = 1 # number of processes
# Output folder to store mesh related files
OUTPUT_FILES = OUTPUT_FILES/attenuation_on
#-----------------------------------------------------------
#
# Mesh
#
#-----------------------------------------------------------
# Partitioning algorithm for decompose_mesh
PARTITIONING_TYPE = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
# number of control nodes per element (4 or 9)
NGNOD = 9
# location to store the mesh
database_filename = OUTPUT_FILES/attenuation_on/database.bin
#-----------------------------------------------------------
#
# Receivers
#
#-----------------------------------------------------------
# use an existing STATION file found in ./DATA or create a new one from the receiver positions below in this Par_file
use_existing_STATIONS = .false.
# number of receiver sets (i.e. number of receiver lines to create below)
nreceiversets = 2
# orientation
anglerec = 0.d0 # angle to rotate components at receivers
rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
# first receiver set (repeat these 6 lines and adjust nreceiversets accordingly)
nrec = 3 # number of receivers
xdeb = 2200. # first receiver x in meters
zdeb = 2200. # first receiver z in meters
xfin = 2800. # last receiver x in meters (ignored if only one receiver)
zfin = 2200. # last receiver z in meters (ignored if only one receiver)
record_at_surface_same_vertical = .true. # receivers inside the medium or at the surface (z values are ignored if this is set to true, they are replaced with the topography height)
# second receiver set
nrec = 3 # number of receivers
xdeb = 2500. # first receiver x in meters
zdeb = 2500. # first receiver z in meters
xfin = 2500. # last receiver x in meters (ignored if only one receiver)
zfin = 1900. # last receiver z in meters (ignored if only one receiver)
record_at_surface_same_vertical = .false. # receivers inside the medium or at the surface (z values are ignored if this is set to true, they are replaced with the topography height)
# filename to store stations file
stations_filename = OUTPUT_FILES/attenuation_on/STATIONS
#-----------------------------------------------------------
#
# Velocity and density models
#
#-----------------------------------------------------------
# number of model materials
nbmodels = 1
# available material types (see user manual for more information)
# acoustic: model_number 1 rho Vp 0 0 0 QKappa 9999 0 0 0 0 0 0 (for QKappa use 9999 to ignore it)
# elastic: model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0 (for QKappa and Qmu use 9999 to ignore them)
# anisotropic: model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 0 QKappa Qmu
# anisotropic in AXISYM: model_number 2 rho c11 c13 c15 c33 c35 c55 c12 c23 c25 c22 QKappa Qmu
# poroelastic: model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu
# tomo: model_number -1 0 0 A 0 0 0 0 0 0 0 0 0 0
#
# note: When viscoelasticity or viscoacousticity is turned on,
# the Vp and Vs values that are read here are the UNRELAXED ones i.e. the values at infinite frequency
# unless the READ_VELOCITIES_AT_f0 parameter above is set to true, in which case they are the values at frequency f0.
#
# Please also note that Qmu is always equal to Qs, but Qkappa is in general not equal to Qp.
# To convert one to the other see doc/Qkappa_Qmu_versus_Qp_Qs_relationship_in_2D_plane_strain.pdf and
# utils/attenuation/conversion_from_Qkappa_Qmu_to_Qp_Qs_from_Dahlen_Tromp_959_960.f90.
1 1 2700.d0 3000.d0 1732.051d0 0 0 100 50 0 0 0 0 0 0
# external tomography file
TOMOGRAPHY_FILE = ./DATA/tomo_file.xyz
# use an external mesh created by an external meshing tool or use the internal mesher
read_external_mesh = .false.
#-----------------------------------------------------------
#
# PARAMETERS FOR EXTERNAL MESHING
#
#-----------------------------------------------------------
# data concerning mesh, when generated using third-party app (more info in README)
# (see also absorbing_conditions above)
mesh_file = ./DATA/mesh_file # file containing the mesh
nodes_coords_file = ./DATA/nodes_coords_file # file containing the nodes coordinates
materials_file = ./DATA/materials_file # file containing the material number for each element
free_surface_file = ./DATA/free_surface_file # file containing the free surface
axial_elements_file = ./DATA/axial_elements_file # file containing the axial elements if AXISYM is true
absorbing_surface_file = ./DATA/absorbing_surface_file # file containing the absorbing surface
acoustic_forcing_surface_file = ./DATA/MSH/Surf_acforcing_Bottom_enforcing_mesh # file containing the acoustic forcing surface
absorbing_cpml_file = ./DATA/absorbing_cpml_file # file containing the CPML element numbers
tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
#-----------------------------------------------------------
#
# PARAMETERS FOR INTERNAL MESHING
#
#-----------------------------------------------------------
# file containing interfaces for internal mesh
interfacesfile = topography.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
xmax = 4000.d0 # abscissa of right side of the model
nx = 80 # number of elements along X
STACEY_ABSORBING_CONDITIONS = .false.
# absorbing boundary parameters (see absorbing_conditions above)
absorbbottom = .true.
absorbright = .true.
absorbtop = .true.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral-element mesh
nbregions = 1 # then set below the different regions and model number for each region
# format of each line: nxmin nxmax nzmin nzmax material_number
1 80 1 60 1
#-----------------------------------------------------------
#
# DISPLAY PARAMETERS
#
#-----------------------------------------------------------
# meshing output
output_grid_Gnuplot = .false. # generate a GNUPLOT file containing the grid, and a script to plot it
output_grid_ASCII = .false. # dump the grid in an ASCII text file consisting of a set of X,Y,Z points or not
The only difference is the material line: QKappa=100 Qmu=50 — these finite quality
factors activate viscoelastic attenuation. Lower Q values mean stronger attenuation, making
the effect clearly visible in the seismograms.
Note
Qmu is always equal to Qs. QKappa is in general not equal to Qp; see
doc/Qkappa_Qmu_versus_Qp_Qs_relationship_in_2D_plane_strain.pdf for the conversion.
Topography file¶
topography.dat¶
#
# number of interfaces
#
2
#
# for each interface below, we give the number of points and then x,z for each point
#
#
# interface number 1 (bottom of the mesh)
#
2
0 0
5000 0
# interface number 2 (topography, top of the mesh)
#
2
0 3000
5000 3000
#
# for each layer, we give the number of spectral elements in the vertical direction
#
#
# layer number 1 (bottom layer)
#
60
The topography file defines the mesh interfaces. The model has a single layer spanning a 4000 m × 3000 m domain, discretized into 80 × 60 spectral elements.
Running xmeshfem2D¶
Generate the mesh database for each variant:
xmeshfem2D -p Par_File_attenuation_off
xmeshfem2D -p Par_File_attenuation_on
Each directory should now contain database.bin and STATIONS files, along with VTK
files for visualization.
Defining sources¶
Both runs share the same source, defined in a single YAML file. For a full description of source parameters refer to Source Description.
source.yaml¶
number-of-sources: 1
sources:
- force:
x : 2500.0
z : 2500.0
source_surf: false
angle : 0.0
vx : 0.0
vz : 0.0
Ricker:
factor: 1e10
tshift: 0.0
f0: 12.5
A single point force source is placed at the center of the domain (2500 m, 2500 m) and excited with a Ricker wavelet at a peak frequency of 12.5 Hz.
Configuring the solver¶
Without attenuation¶
specfem_config_attenuation_off.yaml¶
parameters:
header:
## Header information is used for logging. It is good practice to give your simulations explicit names
title: Isotropic Elastic simulation # name for your simulation
# A detailed description for your simulation
description: |
Material systems : Elastic domain (1)
Interfaces : None
Sources : Force source (1)
Boundary conditions : Neumann BCs on all edges
simulation-setup:
## quadrature setup
quadrature:
quadrature-type: GLL4
## Solver setup
solver:
time-marching:
time-scheme:
type: Newmark
dt: 1.1e-3
nstep: 1600
simulation-mode:
forward:
writer:
seismogram:
format: "ascii"
directory: "OUTPUT_FILES/attenuation_off/results"
receivers:
stations: "OUTPUT_FILES/attenuation_off/STATIONS"
angle: 0.0
seismogram-type:
- velocity
nstep_between_samples: 1
## Runtime setup
run-setup:
number-of-processors: 1
number-of-runs: 1
## databases
databases:
mesh-database: "OUTPUT_FILES/attenuation_off/database.bin"
## sources
sources: "source.yaml"
## log-file configuration
log-file: "OUTPUT_FILES/attenuation_off/output"
This configuration uses the Newmark time scheme with dt=1.1e-3 s for 1600 steps
(total simulation time ~1.76 s). There is no attenuation section, so the solver
treats the medium as purely elastic.
With attenuation¶
specfem_config_attenuation_on.yaml¶
parameters:
header:
## Header information is used for logging. It is good practice to give your simulations explicit names
title: Isotropic Elastic simulation # name for your simulation
# A detailed description for your simulation
description: |
Material systems : Elastic domain (1)
Interfaces : None
Sources : Force source (1)
Boundary conditions : Neumann BCs on all edges
simulation-setup:
## quadrature setup
quadrature:
quadrature-type: GLL4
## Solver setup
solver:
time-marching:
time-scheme:
type: Newmark
dt: 1.1e-3
nstep: 1600
simulation-mode:
forward:
writer:
seismogram:
format: "ascii"
directory: "OUTPUT_FILES/attenuation_on/results"
attenuation:
enabled: true
reference-frequency: 30.0 Hz
attenuation-frequency-band:
- 3.0 Hz
- 300.0 Hz
receivers:
stations: "OUTPUT_FILES/attenuation_on/STATIONS"
angle: 0.0
seismogram-type:
- velocity
nstep_between_samples: 1
## Runtime setup
run-setup:
number-of-processors: 1
number-of-runs: 1
## databases
databases:
mesh-database: "OUTPUT_FILES/attenuation_on/database.bin"
## sources
sources: "source.yaml"
## log-file configuration
log-file: "OUTPUT_FILES/attenuation_on/output"
The attenuation-enabled configuration adds the attenuation section (lines 33–38):
enabled: Set to
trueto activate viscoelastic attenuation.reference-frequency: The frequency at which the velocity model is specified. Set to 30 Hz, within the source frequency band.
attenuation-frequency-band: The band
[f_min, f_max]over which the standard linear solid (SLS) mechanism approximates the target Q. Set to[3 Hz, 300 Hz]to bracket the 12.5 Hz source peak. Choose this band to cover the dominant frequencies of your source.
Running the solver¶
Run both variants back-to-back (or in parallel if resources allow):
specfem 2d -p specfem_config_attenuation_off.yaml
specfem 2d -p specfem_config_attenuation_on.yaml
Note
Make sure either you are in the executable directory of SPECFEM++ or the
executable directory is added to your PATH.
Each run writes seismograms to its respective OUTPUT_FILES/<variant>/results/
directory.
Visualizing the geometry¶
The following script plots the model domain, source location, and receiver positions to verify the setup before running the solver.
"""Plot the benchmark geometry: domain box, source, and station locations."""
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import yaml
TOPOGRAPHY_FILE = "topography.dat"
PAR_FILE = "Par_File_attenuation_on"
SOURCE_FILE = "source.yaml"
STATIONS_FILE = "OUTPUT_FILES/attenuation_on/STATIONS"
OUTPUT_FILE = "OUTPUT_FILES/results/geometry.png"
# ---------------------------------------------------------------------------
# Parse Par_File for xmin/xmax
# ---------------------------------------------------------------------------
def parse_par_file(path):
params = {}
with open(path) as f:
for line in f:
line = line.split("#")[0].strip()
if "=" in line:
key, _, val = line.partition("=")
params[key.strip()] = val.strip().split()[0]
xmin = float(params.get("xmin", "0").rstrip("d0").rstrip("d"))
xmax = float(params.get("xmax", "5000").rstrip("d0").rstrip("d"))
return xmin, xmax
# ---------------------------------------------------------------------------
# Parse topography.dat — returns list of interfaces, each a list of (x,z)
# ---------------------------------------------------------------------------
def parse_topography(path):
with open(path) as f:
lines = [l1.split("#")[0].strip() for l1 in f]
lines = [l2 for l2 in lines if l2]
idx = 0
n_interfaces = int(lines[idx])
idx += 1
interfaces = []
for _ in range(n_interfaces):
n_pts = int(lines[idx])
idx += 1
pts = []
for _ in range(n_pts):
x, z = map(float, lines[idx].split())
idx += 1
pts.append((x, z))
interfaces.append(pts)
return interfaces
# ---------------------------------------------------------------------------
# Parse source.yaml — returns list of (x, z) source positions
# ---------------------------------------------------------------------------
def parse_sources(path):
with open(path) as f:
doc = yaml.safe_load(f)
sources = []
for src in doc.get("sources", []):
for src_type, params in src.items():
x = float(params.get("x", 0))
z = float(params.get("z", 0))
sources.append((x, z))
return sources
# ---------------------------------------------------------------------------
# Parse STATIONS file — columns: name network x z ...
# ---------------------------------------------------------------------------
def parse_stations(path):
stations = []
with open(path) as f:
for line in f:
line = line.strip()
if not line or line.startswith("#"):
continue
parts = line.split()
name = parts[0]
network = parts[1]
x = float(parts[2])
z = float(parts[3])
stations.append((name, network, x, z))
return stations
# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------
def main():
xmin, xmax = parse_par_file(PAR_FILE)
interfaces = parse_topography(TOPOGRAPHY_FILE)
sources = parse_sources(SOURCE_FILE)
stations = parse_stations(STATIONS_FILE)
bottom_pts = interfaces[0]
top_pts = interfaces[-1]
def interp_z(pts, x):
xs = [p[0] for p in pts]
zs = [p[1] for p in pts]
return float(np.interp(x, xs, zs))
bot_xs = np.array([p[0] for p in bottom_pts])
bot_zs = np.array([p[1] for p in bottom_pts])
top_xs = np.array([p[0] for p in top_pts])
top_zs = np.array([p[1] for p in top_pts])
mask_bot = (bot_xs >= xmin) & (bot_xs <= xmax)
mask_top = (top_xs >= xmin) & (top_xs <= xmax)
bot_x = np.concatenate([[xmin], bot_xs[mask_bot], [xmax]])
bot_z = np.concatenate(
[[interp_z(bottom_pts, xmin)], bot_zs[mask_bot], [interp_z(bottom_pts, xmax)]]
)
top_x = np.concatenate([[xmin], top_xs[mask_top], [xmax]])
top_z = np.concatenate(
[[interp_z(top_pts, xmin)], top_zs[mask_top], [interp_z(top_pts, xmax)]]
)
domain_x = np.concatenate([bot_x, top_x[::-1], [bot_x[0]]])
domain_z = np.concatenate([bot_z, top_z[::-1], [bot_z[0]]])
fig, ax = plt.subplots(figsize=(8, 5))
ax.fill(domain_x, domain_z, color="#e8f4f8", zorder=0)
ax.plot(domain_x, domain_z, color="steelblue", linewidth=1.5, zorder=1)
ax.plot(top_x, top_z, color="sienna", linewidth=2, label="Free surface", zorder=2)
side_z = [interp_z(bottom_pts, xmin), interp_z(top_pts, xmin)]
ax.plot(
[xmin, xmin],
side_z,
color="gray",
linewidth=1.5,
linestyle="--",
label="Absorbing boundary",
zorder=2,
)
ax.plot(
[xmax, xmax],
[interp_z(bottom_pts, xmax), interp_z(top_pts, xmax)],
color="gray",
linewidth=1.5,
linestyle="--",
zorder=2,
)
ax.plot(bot_x, bot_z, color="gray", linewidth=1.5, linestyle="--", zorder=2)
sx = [s[2] for s in stations]
sz = [s[3] for s in stations]
sc = ax.scatter(
sx, sz, marker="v", color="royalblue", s=80, zorder=5, label="Stations"
)
sc.set_clip_on(False)
z_surface = interp_z(top_pts, xmin)
surface_count = 0
for name, network, x, z in stations:
at_surface = abs(z - z_surface) < 1.0
if at_surface:
xytext = (0, 6) if surface_count % 2 == 0 else (0, -12)
surface_count += 1
else:
xytext = (4, -12)
ax.annotate(
f"{network}.{name}",
(x, z),
textcoords="offset points",
xytext=xytext,
fontsize=6,
color="royalblue",
ha="center",
annotation_clip=False,
)
for x, z in sources:
ax.scatter(x, z, marker="*", color="crimson", s=200, zorder=6, label="Source")
handles, labels = ax.get_legend_handles_labels()
by_label = dict(zip(labels, handles))
ax.legend(
by_label.values(),
by_label.keys(),
loc="lower left",
fontsize=8,
fancybox=False,
framealpha=1.0,
)
ax.set_xlabel("x (m)")
ax.set_ylabel("z (m)")
ax.set_aspect("equal")
ax.set_xlim(xmin, xmax)
ax.set_ylim(min(bot_z), max(top_z))
for spine in ax.spines.values():
spine.set_visible(False)
fig.tight_layout()
fig.savefig(OUTPUT_FILE, dpi=150)
plt.close(fig)
print(f"Saved geometry plot to {OUTPUT_FILE}")
if __name__ == "__main__":
main()
python plot_geometry.py
Expected geometry:
Model geometry showing the 4000 m × 3000 m domain, absorbing boundaries (dashed gray), the source (red star), and receiver locations (blue triangles).¶
Visualizing seismograms¶
The following Python script reads both sets of seismograms with obspy and produces
comparison plots for the X and Z components.
import glob
import os
import numpy as np
import obspy
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
def get_traces(directory):
traces = []
files = glob.glob(directory + "/*.sem*")
for filename in files:
station_name = os.path.splitext(filename)[0]
network, station, location, channel = station_name.split("/")[-1].split(".")
trace = np.loadtxt(filename, delimiter=" ")
starttime = trace[0, 0]
dt = trace[1, 0] - trace[0, 0]
traces.append(
obspy.Trace(
trace[:, 1],
{
"network": network,
"station": station,
"location": location,
"channel": channel,
"starttime": starttime,
"delta": dt,
},
)
)
return obspy.Stream(traces)
def plot_comparison(stream, stream_att, component, outfile):
sel = stream.select(component=component)
sel_att = stream_att.select(component=component)
stations = sorted(set(tr.stats.station for tr in sel))
fig, axes = plt.subplots(
len(stations), 1, figsize=(10, 2 * len(stations)), sharex=True
)
if len(stations) == 1:
axes = [axes]
for ax, station in zip(axes, stations):
tr = sel.select(station=station)[0]
tr_att = sel_att.select(station=station)[0]
ax.plot(tr.times(), tr.data, label="attenuation off", color="black")
ax.plot(
tr_att.times(),
tr_att.data,
label="attenuation on",
color="red",
linestyle="--",
)
ax.set_ylabel(station)
ax.legend(loc="upper right", fontsize=6)
axes[-1].set_xlabel("Time (s)")
fig.suptitle(f"Component {component}: attenuation on vs. off")
fig.tight_layout()
fig.savefig(outfile, dpi=150)
plt.close(fig)
stream = get_traces("OUTPUT_FILES/attenuation_off/results")
stream_att = get_traces("OUTPUT_FILES/attenuation_on/results")
os.makedirs("OUTPUT_FILES/results", exist_ok=True)
plot_comparison(stream, stream_att, "X", "OUTPUT_FILES/results/traces_X.png")
plot_comparison(stream, stream_att, "Z", "OUTPUT_FILES/results/traces_Z.png")
python plot_traces.py
The comparison plots are saved to OUTPUT_FILES/results/.
Expected results¶
X-component (horizontal) velocity seismograms comparing the elastic run (black) with the viscoelastic run (red dashed).¶
Z-component (vertical) velocity seismograms comparing the elastic run (black) with the viscoelastic run (red dashed).¶
In both figures you should observe:
Amplitude decay: the attenuating run (red) has lower peak amplitudes, especially at later arrival times and higher-frequency content.
Phase shift / velocity dispersion: attenuation is physically linked to dispersion (Kramers–Kronig relations), so the attenuating wavefield arrives very slightly earlier or later depending on frequency.
Waveform broadening: high-frequency energy is preferentially absorbed, causing the waveform to appear smoother in the attenuating run.
About attenuation in SPECFEM++¶
SPECFEM++ models viscoelastic attenuation using the standard linear solid (SLS)
mechanism, implemented via memory variables (also called internal variables). The
approach closely follows Komatitsch & Tromp (1999). The quality factors QKappa
and Qmu in the Par_File control compressional and shear attenuation respectively.
Setting either to 9999 effectively disables that component of attenuation.
The attenuation-frequency-band in the solver configuration determines over which
bandwidth the SLS approximation is optimized. For best accuracy, set this band to
encompass the dominant frequency content of your source.
For full documentation of attenuation parameters refer to SPECFEM++ Parameter Documentation.