Meshes, basis and C-PML¶
PropSEM builds all of this from shape, n_order, topography and abcn.
The builders are public for driving SEMPropagator directly:
a mesh of your own, a C-PML profile tuned by hand, or an equation PropSEM
does not drive.
GLL basis¶
sweep.sem.gll_nodes_weights ¶
gll_nodes_weights(
n_order: int, *, dtype: dtype = torch.float32, device: device | str = "cpu"
) -> tuple[torch.Tensor, torch.Tensor]
N+1 GLL nodes and weights on [-1, 1] as torch tensors.
Always computed in NumPy float64 then cast — sub-ULP precision matters
for stiffness assembly. n_order is the polynomial degree N
(SPECFEM default N=4 → 5 nodes).
sweep.sem.lagrange_derivative_matrix ¶
Derivative matrix D[i, j] = ℓ'_j(x_i) for Lagrange basis at nodes.
for distinct nodes {x_k},
ℓ'j(x_i) = (1 / (x_i - x_j)) · Π (x_i - x_k) / (x_j - x_k) (i ≠ j) ℓ'j(x_j) = Σ 1 / (x_j - x_k)
For GLL nodes there's an even simpler closed form via P_N values (see Canuto Appendix B), but we use the general one for clarity.
2-D meshes¶
Quad4Mesh is the flat, axis-aligned mesh with a constant Jacobian; Quad9Mesh
is the curved one (biquadratic elements) that follows a topography.
sweep.sem.Quad4Mesh ¶
Quad4Mesh(
n_order: int,
nex: int,
nez: int,
x_max: float,
z_max: float,
xi: Tensor,
w_ref: Tensor,
deriv: Tensor,
dx_e: float,
dz_e: float,
jacobian_det: float,
ibool: Tensor,
x_global: Tensor,
z_global: Tensor,
n_global: int,
x_e: Tensor,
z_e: Tensor,
)
Uniform axis-aligned rectangular mesh with GLL node connectivity.
sweep.sem.build_quad4_mesh ¶
build_quad4_mesh(
nex: int,
nez: int,
x_max: float,
z_max: float,
n_order: int = 4,
*,
dtype: dtype = torch.float32,
device: device | str = "cpu"
) -> Quad4Mesh
Build a uniform Cartesian quad mesh with GLL connectivity.
Parameters¶
nex, nez : number of elements in x, z. x_max, z_max : physical domain extents (z increases downward, top is z=0). n_order : polynomial order N. Default 4 → 5×5 GLL nodes per element.
sweep.sem.Quad9Mesh ¶
Quad9Mesh(
n_order: int,
nex: int,
nez: int,
x_max: float,
z_max: float,
xi: Tensor,
w_ref: Tensor,
deriv: Tensor,
control_xy: Tensor,
inv_J: Tensor,
det_J: Tensor,
w_jac: Tensor,
ibool: Tensor,
x_global: Tensor,
z_global: Tensor,
n_global: int,
x_e: Tensor,
z_e: Tensor,
)
Deformed-quad SEM mesh with Q2 geometry and per-GLL Jacobian.
sweep.sem.build_axis_aligned_quad9 ¶
build_axis_aligned_quad9(
nex: int,
nez: int,
x_max: float,
z_max: float,
*,
n_order: int = 4,
dtype: dtype = torch.float32,
device: device | str = "cpu"
) -> Quad9Mesh
Affine Quad9 mesh = Quad4 with Q2 control nodes; used for degeneracy tests.
sweep.sem.build_topography_quad9 ¶
build_topography_quad9(
nex: int,
nez: int,
x_max: float,
z_bottom: float,
z_top_fn: Callable[[Tensor], Tensor],
*,
n_order: int = 4,
dtype: dtype = torch.float32,
device: device | str = "cpu"
) -> Quad9Mesh
Build a conforming Quad9 mesh with the top edge following z_top_fn(x).
The bottom edge is flat at z = z_bottom. Each column of elements
is graded linearly in z between the (curved) top and the (flat)
bottom — the simplest conforming mesh that exactly carries the
topography along the free surface.
z_top_fn must accept a 1-D torch tensor of x coordinates and
return the surface z (in the same SEM down-positive convention,
so z_top_fn(x) ≥ 0 everywhere and z_top_fn(x) < z_bottom).
The element-corner / edge-midpoint x-positions are uniformly spaced in x; z-positions inside each column are linearly graded.
sweep.sem.quad9_from_control_nodes ¶
quad9_from_control_nodes(
control_xy: Tensor,
*,
n_order: int = 4,
dtype: dtype = torch.float32,
device: device | str = "cpu"
) -> Quad9Mesh
Build a deformed mesh from a per-element control-node tensor.
control_xy[ez, ex, i, j, :] is the (x, z) coordinate of the
(i, j)-th Q2 control node of element (ez, ex). Shared GLL nodes on
inter-element edges are de-duplicated by the analytic ibool indexing and
take their coordinate from their lowest-numbered element copy (the mesh is
assumed conforming, so the copies agree up to rounding); the build checks
only det_J > 0 (no inversion), then assembles Jacobian + ibool + global
node coords.
sweep.sem.boundary_element_edges ¶
Return element indices whose edges touch each physical boundary.
'left' (x=0), 'right' (x=x_max),
'top' (z=0), 'bottom' (z=z_max).
Values are 1-D int64 tensors of element-flat indices e = ez*nex + ex.
3-D meshes¶
Hex8Mesh and Hex27Mesh are the 3-D counterparts of Quad4Mesh and
Quad9Mesh.
sweep.sem.Hex8Mesh ¶
Hex8Mesh(
n_order: int,
nex: int,
ney: int,
nez: int,
x_max: float,
y_max: float,
z_max: float,
xi: Tensor,
w_ref: Tensor,
deriv: Tensor,
dx_e: float,
dy_e: float,
dz_e: float,
jacobian_det: float,
ibool: Tensor,
x_global: Tensor,
y_global: Tensor,
z_global: Tensor,
n_global: int,
x_e: Tensor,
y_e: Tensor,
z_e: Tensor,
layout: Hex8Layout | None = None,
)
Uniform axis-aligned rectangular hex mesh with GLL node connectivity.
sweep.sem.build_hex8_mesh ¶
build_hex8_mesh(
nex: int,
ney: int,
nez: int,
x_max: float,
y_max: float,
z_max: float,
n_order: int = 4,
*,
dtype: dtype = torch.float32,
device: device | str = "cpu",
geometry: bool = True,
ibool32: bool = False
) -> Hex8Mesh
Build a uniform Cartesian hex mesh with GLL connectivity.
Parameters¶
nex, ney, nez : number of elements in x, y, z.
x_max, y_max, z_max : physical domain extents (z increases downward,
top is z=0).
n_order : polynomial order N. Default 4 → 5×5×5 GLL nodes per element.
geometry : with False the six coordinate arrays (x_e/y_e/z_e,
x/y/z_global -- 35 of the mesh's 50.7 B per node) are left
None and only the connectivity is built. For a consumer that
reads ibool and the element counts alone, such as a
domain-decomposition strip whose coordinates come from the global
:class:Hex8Layout.
ibool32 : store the connectivity as int32 (the width the kernels take).
Opt-in, and only for such a consumer: a DD strip otherwise holds an
int64 mesh copy AND the int32 kernel copy, measured 15.5 B per global
node per rank. The single-domain path and the boundary helpers read
Hex8Mesh.ibool expecting the default int64.
sweep.sem.Hex27Mesh ¶
Hex27Mesh(
n_order: int,
nex: int,
ney: int,
nez: int,
x_max: float,
y_max: float,
z_max: float,
xi: Tensor,
w_ref: Tensor,
deriv: Tensor,
control_xyz: Tensor,
inv_J: Tensor,
det_J: Tensor,
w_jac: Tensor,
ibool: Tensor,
x_global: Tensor,
y_global: Tensor,
z_global: Tensor,
n_global: int,
x_e: Tensor,
y_e: Tensor,
z_e: Tensor,
)
Deformed-hex SEM mesh with Q2 geometry and per-GLL 3×3 Jacobian.
sweep.sem.build_axis_aligned_hex27 ¶
build_axis_aligned_hex27(
nex, ney, nez, x_max, y_max, z_max, *, n_order=4, dtype=torch.float32, device="cpu"
) -> Hex27Mesh
Affine Hex27 mesh = Hex8 with Q2 control nodes (degeneracy tests).
sweep.sem.build_topography_hex27 ¶
build_topography_hex27(
nex,
ney,
nez,
x_max,
y_max,
z_bottom,
z_top_fn: Callable[[Tensor, Tensor], Tensor],
*,
n_order=4,
dtype=torch.float32,
device="cpu"
) -> Hex27Mesh
Conforming Hex27 mesh whose top face follows z_top_fn(x, y).
The bottom face is flat at z = z_bottom; each vertical column is graded
linearly in z between the curved top surface and the flat bottom. z_top_fn
takes two equal-shaped tensors of (x, y) and returns the surface z (SEM
down-positive: 0 ≤ z_top < z_bottom). x and y control positions are
uniform.
sweep.sem.hex27_from_control_nodes ¶
hex27_from_control_nodes(
control_xyz: Tensor,
nex: int,
ney: int,
nez: int,
*,
n_order: int = 4,
dtype: dtype = torch.float32,
device: device | str = "cpu"
) -> Hex27Mesh
Assemble Jacobian + ibool + global node coords from element control nodes.
sweep.sem.boundary_element_faces ¶
Element indices whose faces touch each physical boundary.
'left' (x=0), 'right' (x=x_max),
'front' (y=0), 'back' (y=y_max),
'top' (z=0), 'bottom' (z=z_max).
Values are 1-D int64 tensors of element-flat indices
e = (ez*ney + ey)*nex + ex.
Mass matrix¶
sweep.sem.assemble_mass_vector ¶
Assemble the global diagonal mass vector M[g] = Σ_e w_i w_j |J_e| · ρc²-inv.
Phase 1: constant ρ acoustic, so user passes 1/(ρ c²) sampled at GLL points. Returns a (n_global,) tensor.
C-PML profiles¶
One profile per equation family; the equation's build takes it as cpml
(the coupled fluid-solid form takes an acoustic and an elastic one).
sweep.sem.CPMLProfile ¶
CPMLProfile(
sides: tuple[str, ...],
pml_thickness: float,
a_x: Tensor,
b_x: Tensor,
dbx_dx: Tensor,
a_z: Tensor,
b_z: Tensor,
dbz_dz: Tensor,
)
Per-global-node CPML recursion coefficients (double-convolution form).
sweep.sem.build_cpml ¶
build_cpml(
mesh: Quad4Mesh,
*,
dt: float,
sides: tuple[str, ...] = ("left", "right", "bottom"),
pml_thickness_elements: int = 8,
c_max: float,
f_dominant: float,
reflection_target: float = 0.001,
poly: float = 2.0
) -> CPMLProfile
Build a Komatitsch-Martin 2007 double-convolution CPML profile.
Parameters¶
sides
Subset of {"left", "right", "top", "bottom"}. Top is usually
omitted (free surface).
pml_thickness_elements
Number of mesh elements in each PML strip; physical thickness =
pml_thickness_elements · dh_e.
c_max, f_dominant
Peak velocity and dominant source frequency, used to derive
σ_max and α_max.
reflection_target
Target normal-incidence reflection coefficient (drives σ_max via
Collino-Tsogka 2001 formula).
poly
Polynomial taper power; 2 is the conventional quadratic profile.
sweep.sem.CPML3DProfile ¶
CPML3DProfile(
sides: tuple[str, ...],
pml_thickness: float,
a_x: Tensor,
b_x: Tensor,
dbx_dx: Tensor,
a_y: Tensor,
b_y: Tensor,
dby_dy: Tensor,
a_z: Tensor,
b_z: Tensor,
dbz_dz: Tensor,
)
Per-global-node CPML recursion coefficients (3-D double-convolution).
sweep.sem.build_cpml3d ¶
build_cpml3d(
mesh: Hex8Mesh,
*,
dt: float,
sides: tuple[str, ...] = ("left", "right", "front", "back", "bottom"),
pml_thickness_elements: int = 8,
c_max: float,
f_dominant: float,
reflection_target: float = 0.001,
poly: float = 2.0
) -> CPML3DProfile
Build a Komatitsch-Martin 2007 double-convolution 3-D CPML profile.
Parameters¶
sides
Subset of {"left", "right", "front", "back", "top", "bottom"}.
Top is usually omitted (free surface).
pml_thickness_elements
Number of mesh elements in each PML strip.
c_max, f_dominant
Peak velocity and dominant source frequency (drive σ_max / α_max).
reflection_target
Target normal-incidence reflection coefficient (Collino-Tsogka 2001).
poly
Polynomial taper power (2 = conventional quadratic).
sweep.sem.ElasticCPMLProfile ¶
ElasticCPMLProfile(
sides: tuple,
pml_thickness: float,
node_in_pml: Tensor,
Sx0: Tensor,
Sx1: Tensor,
Sx2: Tensor,
zx1: tuple,
zx2: tuple,
Sz0: Tensor,
Sz1: Tensor,
Sz2: Tensor,
xz1: tuple,
xz2: tuple,
La1: Tensor,
La2: Tensor,
La3: Tensor,
La4: Tensor,
l1: tuple,
l2: tuple,
outer_edge_idx: Tensor,
mpml_ratio: float = 0.0,
)
Per-element-GLL SPECFEM-style coefficients for the 2nd-order elastic CPML.
Every field has shape (E, ngll, ngll). Outside the PML the strain
coefficients reduce to a no-op (S?0 = 1, S?1 = S?2 = 0) and the
accel coefficients vanish, so the recursions are safe to run everywhere.
The recursive-convolution coefficients are masked to zero outside the PML so
the (unused) memory branches never accumulate.
mpml_ratio
class-attribute
¶
Convert a string or number to a floating point number, if possible.
sweep.sem.build_elastic_cpml ¶
build_elastic_cpml(
mesh: Quad9Mesh,
*,
dt: float,
sides: tuple = ("left", "right", "bottom"),
pml_thickness_elements: int = 8,
c_max: float,
f_dominant: float,
reflection_target: float = 0.001,
poly: float = 2.0,
mpml_ratio: float = 0.0
) -> ElasticCPMLProfile
Compute the full SPECFEM2D elastic-CPML coefficient set at every GLL node.
SPECFEM recipe (K_MAX = 1, no rotation):
d_max = -ln(R)·(p+1)·c / (2 L); α_max = π f₀
d(s) = d_max·(s/L)^p inward; α(s) = α_max·(1 − s/L)
β = α + d (since κ = 1)
then lik (strain stretch, two calls) and l (accel term).
mpml_ratio (p in Meza-Fajardo & Papageorgiou 2008): the Multiaxial
PML cross-damping fraction. Classical (uniaxial) C-PML damps only normal
to each boundary (d_x in the left/right bands, d_z in the bottom band),
which is weakly reflective/unstable for grazing and surface (Rayleigh) waves in
elastic media. mpml_ratio > 0 adds a small fraction of the normal damping
in the tangential direction — d_x += p·d_z⁰, d_z += p·d_x⁰ — so every
band becomes multiaxial (a "corner" in the region sense). This stabilises the
surface-wave absorption at the cost of the layer no longer being perfectly
matched (a small reflection). Typical p ∈ [0.1, 0.3]; 0 = classical
C-PML (default, exact SPECFEM2D behaviour).
sweep.sem.SimpleElasticCPMLProfile ¶
SimpleElasticCPMLProfile(
sides: tuple,
pml_thickness: float,
a_x: Tensor,
b_x: Tensor,
a_z: Tensor,
b_z: Tensor,
)
FD-style 1-stage CPML profile, per element-GLL node.
Storage: each tensor has shape (E, ngll, ngll). Outside the PML
both a_x and b_x (or a_z, b_z) are zero, so the
recursion m' = a · m + b · ∂ becomes m' = 0 and ∂̃ = ∂.
sweep.sem.build_simple_elastic_cpml ¶
build_simple_elastic_cpml(
mesh: Quad9Mesh,
*,
dt: float,
sides: tuple = ("left", "right", "bottom"),
pml_thickness_elements: int = 8,
c_max: float,
f_dominant: float,
reflection_target: float = 0.001,
poly: float = 2.0
) -> SimpleElasticCPMLProfile
Build the FD-style CPML profile for mesh.
Parameters mirror :func:sweep.equations.pml.setup_pml:
σ_0 = -(n + 1) · c_max · ln(R) / (2 L) α_max = π · f_0 σ(x) = σ_0 · (depth/L)^n α(x) = α_max · (1 − depth/L)
Then per axis:
a = exp(-(σ + α) · dt) b = σ / (σ + α) · (a − 1)
sweep.sem.SimpleElasticCPML3DProfile ¶
SimpleElasticCPML3DProfile(
sides: tuple,
pml_thickness: float,
a_x: Tensor,
b_x: Tensor,
a_y: Tensor,
b_y: Tensor,
a_z: Tensor,
b_z: Tensor,
)
FD-style 1-stage 3-D CPML profile, per element-GLL node (E,g,g,g).
sweep.sem.build_simple_elastic_cpml_3d ¶
build_simple_elastic_cpml_3d(
mesh: Hex8Mesh,
*,
dt: float,
sides: tuple = ("left", "right", "front", "back", "bottom"),
pml_thickness_elements: int = 8,
c_max: float,
f_dominant: float,
reflection_target: float = 0.001,
poly: float = 2.0
) -> SimpleElasticCPML3DProfile
Build the FD-style 3-D CPML profile for mesh (one (a,b) pair per axis).
σ_0 = -(n + 1) · c_max · ln(R) / (2 L) α_max = π · f_0 σ(s) = σ_0 · (depth/L)^n ; α(s) = α_max · (1 − depth/L)