Example 2: Multi-layer Soil Model with Absorbing Boundaries
Overview
This example demonstrates how to create a larger-scale soil domain with multiple layers and absorbing boundary conditions using Perfectly Matched Layers (PML). The model is designed for wave propagation analysis with domain reduction method (DRM) loading.
The model created in this example consists of:
Six distinct soil layers with depth-varying properties
3D brick elements with elastic-isotropic materials
Rayleigh damping applied to control dynamic response
PML absorbing boundaries to prevent wave reflections
DRM loading applied from an external file

Conceptual visualization of the multi-layer soil model with absorbing boundaries
Model Description
Soil Domain:
128m × 128m domain (-64m to 64m in both X and Y directions)
Total depth of 48m
Six distinct soil layers, each 8m thick
Properties that vary with depth (density, wave velocities, and damping)
Materials:
Elastic-isotropic materials with properties calculated from wave velocities
Density (rho) increasing with depth from 2142 to 2162 kg/m³
Shear wave velocity (Vs) increasing with depth from 353 to 599 m/s
P-wave velocity (Vp) increasing with depth from 669 to 1135 m/s
Damping ratio decreasing with depth
Mesh:
Element sizes: 4.0m × 4.0m × 4.0m uniform grid
32 × 32 elements in horizontal directions
2 elements per layer in vertical direction
Damping:
Frequency-dependent Rayleigh damping
Target frequencies: f1=3Hz and f2=15Hz
Damping ratio varying by layer (from 2.1% to 3.0%)
Boundary Conditions:
Perfectly Matched Layer (PML) absorbing boundaries
4 element layers in PML
Rayleigh damping coefficient of 0.95 in PML region
Loading:
Domain Reduction Method (DRM) loading
Input from external H5DRM file
Simulation duration: 30 seconds with 0.01s time step
Domain Reduction Method (DRM) Explained
The Domain Reduction Method (DRM) is a two-step computational approach for simulating seismic wave propagation in a localized region of interest:
Background Analysis: First, a large domain is analyzed to capture the regional seismic wave propagation from the source to the site.
Local Analysis: The results from step 1 are then used as boundary conditions for a smaller domain of interest, allowing for detailed analysis without modeling the entire domain.
Key features and benefits of DRM:
Computational Efficiency: Enables detailed analysis of a small region without the computational cost of modeling the entire wave propagation path from source to site
Wave Input Flexibility: Can input waves from any direction and source mechanism
Realistic Boundary Conditions: Accounts for complex wave propagation effects more accurately than direct force inputs
Localized Refinement: Allows for mesh refinement in the region of interest without needing to refine the entire domain
In this example, DRM is implemented through:
h5pattern = fm.pattern.create_pattern('h5drm',
filepath='drmload.h5drm',
factor=1.0,
crd_scale=1.0,
distance_tolerance=0.01,
do_coordinate_transformation=1,
transform_matrix=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
origin=[0.0, 0.0, 0.0])
fm.drm.set_pattern(h5pattern)
The h5drm
file contains pre-computed wave field information that is applied to the DRM boundary elements, effectively injecting the desired seismic waves into the computational domain.
Perfectly Matched Layer (PML) Explained
PML is an advanced absorbing boundary technique developed to prevent artificial reflections of waves at computational domain boundaries:
fm.drm.addAbsorbingLayer(numLayers=4,
numPartitions=8,
partitionAlgo="kd-tree",
geometry="Rectangular",
rayleighDamping=0.95,
matchDamping=False,
type="PML")
Key features of PML:
Wave Absorption: PML gradually attenuates waves as they enter the boundary region, preventing them from reflecting back into the domain
Angle Independence: Designed to have minimal reflections regardless of wave incidence angle or frequency content
Computational Efficiency: More efficient than extending the physical domain to prevent reflections
Implementation Parameters: -
numLayers
: Controls thickness of the absorbing region (4 elements in this example) -rayleighDamping
: Sets the damping coefficient in the PML region (0.95 in this example)
PML works by creating an artificial absorbing material whose damping properties gradually increase with distance from the computational domain boundary. This gradual transition ensures waves are absorbed without creating spurious reflections at the PML interface.
Absorbing Layer Parameters in Detail
The addAbsorbingLayer
method creates boundary regions that absorb outgoing waves to prevent them from reflecting back into the computational domain. Let’s examine each parameter in detail:
numLayers (int):
Defines the thickness of the absorbing boundary in terms of element layers
More layers typically provide better absorption but increase computational cost
Recommended values range from 3-8 layers, with 4 being a good balance in many cases
In this example, 4 layers effectively absorb most seismic waves while maintaining efficiency
numPartitions (int):
Controls how the boundary regions are divided for parallel processing
Higher values can improve load balancing but may increase overhead
Value should typically scale with the number of available CPU cores
In this example, 8 partitions are used for efficient parallel computation
partitionAlgo (str):
Specifies the algorithm used to divide the model for parallel processing
Currently implemented options:
“kd-tree”: Spatially balanced partitioning (used in this example)
“metis”: Graph-based partitioning for complex geometries
“kd-tree” is generally recommended for soil models as it provides good load balancing
geometry (str):
Defines the shape of the computational domain and absorbing boundaries
Currently only “Rectangular” is fully implemented
Future implementations will include:
“Cylindrical”: For cylindrical domains
“Spherical”: For spherical domains
The “Rectangular” option is appropriate for this layered soil model with a rectangular footprint
rayleighDamping (float):
Controls the strength of the Rayleigh damping in the absorbing layer
Ranges from 0.0 (no damping) to 1.0 (maximum damping)
Higher values (like 0.95 in this example) absorb more energy but may cause numerical issues
This parameter sets a scalar multiplier for the damping coefficients that increases with distance into the boundary
The damping strength gradually increases from the inner boundary (computational domain) to the outer boundary
matchDamping (bool):
When True, matches the damping properties between the computational domain and absorbing layers
When False (as in this example), uses the damping properties defined specifically for the absorbing layer
Setting to False allows for stronger damping in the absorbing layers than in the computational domain
This parameter is particularly relevant when the computational domain already has some level of damping
type (str):
Specifies the type of absorbing boundary to use
Currently implemented options:
“PML”: Perfectly Matched Layer (used in this example)
“Rayleigh”: Standard Rayleigh damping layer
“ASDA”: Combined PML and Rayleigh approach
Comparing Absorbing Layer Types:
PML (Perfectly Matched Layer):
More effective at absorbing waves at any angle of incidence
Better performance for complex wave fields
More computationally intensive
Preferred for high-precision wave propagation simulations
Rayleigh damping layers:
Simpler implementation
Less computationally demanding
More numerically stable than PML
May require more layers to achieve the same absorption quality
Recommended for models with stability concerns
ASDA (Combined approach):
Combines benefits of both PML and Rayleigh approaches
Can provide better stability than pure PML
More computationally intensive than Rayleigh alone
In this example, PML is chosen because it provides superior absorption characteristics for seismic wave propagation problems, particularly when waves may approach the boundary at various angles. However, if you experience numerical instabilities, switching to “Rayleigh” type may help stabilize the solution at the cost of some absorption efficiency.
Frequency Rayleigh Damping in Femora
In this example, each soil layer is assigned its own frequency-dependent Rayleigh damping:
damp = fm.damping.create_damping("frequency rayleigh", dampingFactor=xi_s, f1=3, f2=15)
Frequency Rayleigh damping is a more intuitive approach to specifying damping than conventional Rayleigh damping. Instead of directly setting the mass and stiffness proportional coefficients (αₘ and βₖ), users can specify:
dampingFactor: The target damping ratio ξ (between 0 and 1) - in our example, this varies by layer (xi_s)
f1: Lower bound target frequency in Hz (3 Hz in this model)
f2: Upper bound target frequency in Hz (15 Hz in this model)
How Frequency Rayleigh Damping Works:
Rayleigh damping combines mass-proportional and stiffness-proportional damping through the equation:
Where: - [C] is the damping matrix - [M] is the mass matrix - [K] is the stiffness matrix - αₘ and βₖ are the mass and stiffness proportional coefficients
Traditional Rayleigh damping requires specifying αₘ and βₖ directly, which is not intuitive. In frequency Rayleigh damping, Femora automatically calculates these coefficients from more practical inputs:
This approach is implemented in Femora because:
It’s more intuitive for geotechnical engineers who think in terms of frequency ranges and damping ratios
It automatically ensures the damping ratio equals the target value exactly at the two specified frequencies
It simplifies the process of creating realistic energy dissipation in dynamic models
It allows for frequency-calibrated damping that properly accounts for the predominant frequencies in seismic analysis
Mathematical Background:
The damping ratio at any frequency ω is given by:
By selecting two target frequencies (ω₁ and ω₂) and setting the damping ratio to the desired value (ξ) at these frequencies, we can solve for the Rayleigh coefficients αₘ and βₖ. This calibration ensures that the damping ratio is exactly ξ at frequencies f₁ and f₂, and varies gradually for other frequencies.
Region Assignment in Layered Soil Model
In this example, each soil layer is assigned to its own region with specific damping properties:
# Create frequency-dependent damping for this layer
damp = fm.damping.create_damping("frequency rayleigh", dampingFactor=xi_s, f1=3, f2=15)
# Create region with this layer's specific damping
reg = fm.region.create_region("elementRegion", damping=damp)
# Assign the region when creating the mesh part for this layer
fm.meshPart.create_mesh_part("Volume mesh", "Uniform Rectangular Grid",
user_name=f"Layer{layer['layer']}",
element=ele,
region=reg,
# ... other parameters ...)
Region Assignment Process:
For each layer in the soil profile, we:
Create a layer-specific damping object with damping ratio corresponding to that layer
Create a new element region and associate it with this damping
Assign this region to the mesh part when creating the layer geometry
This workflow allows for:
Layer-specific damping: Each soil layer can have different damping properties based on its material characteristics
Depth-dependent behavior: Properties can vary with depth, reflecting real soil conditions
Independent energy dissipation: Each region’s response to dynamic loading is calibrated separately
Global Region and Region Management in Femora
All Femora models have a “global region” that serves as the default assignment when no explicit region is provided.
Region System Architecture in Femora:
Femora implements a sophisticated region management system with these key components:
Global Region:
Automatically created when Femora initializes
Has tag 0 and name “GlobalRegion”
Serves as a fallback when no specific region is assigned
Can have its own damping properties that apply model-wide
Element Regions:
Created explicitly with
fm.region.create_region("elementRegion", ...)
Can have specific damping assigned
Can be associated with specific elements or element groups
Region Assignment Rules:
When creating a mesh part, a region must be specified or the global region is used
The statement
region=reg
explicitly assigns a specific region to that mesh partWithout this assignment, all elements would use the global region’s properties
Implementation Details:
The global region is implemented through a singleton pattern to ensure there is always exactly one global region with tag 0:
def __new__(cls, tag=None, damping: Type[DampingBase] = None):
if tag in cls._regions:
return cls._regions[tag]
if tag == 0:
cls._global_region = super().__new__(cls)
cls._regions[0] = cls._global_region
return cls._global_region
When creating a mesh part without specifying a region:
# This would implicitly use the global region
fm.meshPart.create_mesh_part("Volume mesh", "Uniform Rectangular Grid",
user_name=f"SomeLayer",
element=ele)
In our soil model example, explicit region creation ensures that each soil layer has its own properly configured damping behavior rather than using global defaults. If we had omitted the region assignment:
# Without explicit region, would use global region
fm.meshPart.create_mesh_part("Volume mesh", "Uniform Rectangular Grid",
user_name=f"Layer{layer['layer']}",
element=ele)
Then all soil layers would have inherited the same damping properties from the global region, which would not accurately represent the depth-varying damping characteristics of real soil profiles.
This region-based approach allows Femora to implement element-level, region-level, or global-level damping in a consistent way, giving users flexibility to model complex systems with varying properties in different zones while maintaining a simple default behavior through the global region.
Step-by-Step Code Walkthrough
Let’s examine the key components of the Example 2 code:
Setting up the Environment and Layer Properties:
import femora as fm import os # Change directory to the current file's location os.chdir(os.path.dirname(os.path.abspath(__file__))) # Layer properties as a dictionary layers = [ {"layer": 1, "rho": 2142.0500, "vp": 669.0500, "vs": 353.1000, "xi_s": 0.0296, "xi_p": 0.0148, "thickness": 8}, # ... more layers ]
The code begins by importing the Femora library and defining six soil layers with varying properties:
rho
: Material density (kg/m³)vp
: P-wave velocity (m/s)vs
: S-wave velocity (m/s)xi_s
,xi_p
: Damping ratios for shear and pressure wavesthickness
: Layer thickness (8m each)
Defining Rayleigh Damping Parameters:
pi = 3.1415926 f1 = 15; # Hz f2 = 3; # Hz w1 = 2*pi*f1 w2 = 2*pi*f2
Rayleigh damping requires two target frequencies for calibration. These frequencies (3Hz and 15Hz) cover the primary frequency range of interest for seismic analysis.
Creating the Soil Layers:
index = 0 for layer in layers[::-1]: # Calculate elastic properties from wave velocities vp_vs_ratio = (vp / vs) ** 2 nu = (vp_vs_ratio - 2) / (2 * (vp_vs_ratio - 1)) # Poisson's ratio E = 2 * rho * (vs ** 2) * (1 + nu) # Young's modulus # Create material, element, and damping objects mat = fm.material.create_material("nDMaterial", "ElasticIsotropic", user_name=f"Layer{layer['layer']}", E=E, nu=nu, rho=rho) # ... create elements and mesh
This loop processes layers from bottom to top and:
Converts wave velocities to mechanical properties (E and ν)
Creates elastic material models
Applies frequency-dependent Rayleigh damping
Creates a uniform 3D mesh grid for each layer
Model Assembly and Boundary Conditions:
# Create a list of mesh parts mesh_parts = [] for layer in layers: mesh_parts.append(f"Layer{layer['layer']}") fm.assembler.create_section(mesh_parts, num_partitions=4) fm.assembler.Assemble() # Add absorbing boundary layers (PML) fm.drm.addAbsorbingLayer(numLayers=4, numPartitions=8, partitionAlgo="kd-tree", geometry="Rectangular", rayleighDamping=0.95, matchDamping=False, type="PML", )
After creating all layers, the code assembles all mesh parts and adds PML absorbing boundaries.
Loading and Analysis Setup:
h5pattern = fm.pattern.create_pattern('h5drm', filepath='drmload.h5drm', factor=1.0, crd_scale=1.0, distance_tolerance=0.01, do_coordinate_transformation=1, transform_matrix=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0], origin=[0.0, 0.0, 0.0]) fm.drm.set_pattern(h5pattern) fm.drm.createDefaultProcess(finalTime=30, dT=0.01) fm.export_to_tcl(filename="model.tcl") fm.gui()
Finally, the code sets up DRM loading from an H5DRM file, configures the time-domain analysis (30s duration with 0.01s steps), exports the model to a TCL file for analysis, and launches the Femora GUI for visualization.
Code Access
The full source code for this example is available in the Femora repository:
import femora as fm
import os
# change the direcotto the current file
os.chdir(os.path.dirname(os.path.abspath(__file__)))
# Layer properties as a dictionary
layers = [
{"layer": 1, "rho": 2142.0500, "vp": 669.0500, "vs": 353.1000, "xi_s": 0.0296, "xi_p": 0.0148, "thickness": 8},
{"layer": 2, "rho": 2146.2000, "vp": 785.2500, "vs": 414.4000, "xi_s": 0.0269, "xi_p": 0.0134, "thickness": 8},
{"layer": 3, "rho": 2150.3500, "vp": 886.0500, "vs": 467.6000, "xi_s": 0.0249, "xi_p": 0.0125, "thickness": 8},
{"layer": 4, "rho": 2154.4500, "vp": 976.4000, "vs": 515.3000, "xi_s": 0.0234, "xi_p": 0.0117, "thickness": 8},
{"layer": 5, "rho": 2158.6000, "vp": 1059.0500, "vs": 558.9500, "xi_s": 0.0221, "xi_p": 0.0111, "thickness": 8},
{"layer": 6, "rho": 2162.7500, "vp": 1135.6500, "vs": 599.4000, "xi_s": 0.0211, "xi_p": 0.0105, "thickness": 8}
]
pi = 3.1415926
f1 = 15; # Hz
f2 = 3; # Hz
w1 = 2*pi*f1
w2 = 2*pi*f2
index = 0
for layer in layers[::-1]:
rho = layer["rho"]
vp = layer["vp"]
vs = layer["vs"]
xi_s = layer["xi_s"]
xi_p = layer["xi_p"]
thickness = layer["thickness"]
vp_vs_ratio = (vp / vs) ** 2
nu = (vp_vs_ratio - 2) / (2 * (vp_vs_ratio - 1))
E = 2 * rho * (vs ** 2) * (1 + nu)
E = E / 1000
rho = rho / 1000
mat = fm.material.create_material("nDMaterial", "ElasticIsotropic",
user_name=f"Layer{layer['layer']}",
E=E, nu=nu, rho=rho)
ele = fm.element.create_element(element_type="stdBrick",
ndof=3,
material=mat,
b1=0,
b2=0,
b3=0)
damp = fm.damping.create_damping("frequency rayleigh", dampingFactor=xi_s, f1=3, f2=15)
# reg = fm.region.create_region("elementRegion", damping=damp)
zmin = -48 + index * thickness
zmax = zmin + thickness
dx = 4.0
dy = 4.0
dz = 4.0
Nx = int(128 / dx)
Ny = int(128 / dy)
Nz = int(thickness / dz)
reg = fm.region.create_region("elementRegion", damping=damp)
fm.meshPart.create_mesh_part("Volume mesh", "Uniform Rectangular Grid",
user_name=f"Layer{layer['layer']}",
element=ele,
region=reg,
**{'X Min': -64, 'X Max': 64,
'Y Min': -64, 'Y Max': 64,
'Z Min': zmin, 'Z Max': zmax,
'Nx Cells': Nx, 'Ny Cells': Ny, 'Nz Cells': Nz})
index += 1
# create a list of mesh parts
mesh_parts = []
for layer in layers:
mesh_parts.append(f"Layer{layer['layer']}")
fm.assembler.create_section(mesh_parts, num_partitions=4)
fm.assembler.Assemble()
# fm.addAbsorbingLayer(numLayers=10,
# numPartitions=8,
# partitionAlgo="kd-tree",
# geometry="Rectangular",
# rayleighDamping=0.95,
# matchDamping=False,
# type="Rayleigh",
# )
fm.drm.addAbsorbingLayer(numLayers=4,
numPartitions=8,
partitionAlgo="kd-tree",
geometry="Rectangular",
rayleighDamping=0.95,
matchDamping=False,
type="PML",
)
h5pattern = fm.pattern.create_pattern( 'h5drm',
filepath='drmload.h5drm',
factor=1.0,
crd_scale=1.0,
distance_tolerance=0.01,
do_coordinate_transformation=1,
transform_matrix=[1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0],
origin=[0.0, 0.0, 0.0])
fm.drm.set_pattern(h5pattern)
fm.drm.createDefaultProcess(finalTime=30, dT=0.01)
fm.export_to_tcl(filename="model.tcl")
# fm.export_to_vtk("mesh.vtk")
fm.gui()
Example directory:
examples/Example2/
Python script:
examples/Example2/femoramodel.py
Data file:
drmload.h5drm
(required for running the example)
Key Concepts Demonstrated
Creating a multi-layer soil model with depth-varying properties
Calculating material properties from wave velocities
Implementing frequency-dependent Rayleigh damping
Setting up PML absorbing boundaries to prevent wave reflections
Applying DRM loading for seismic wave propagation analysis
Organizing a large-scale soil domain for efficient computation
Assigning regions with specific damping properties
Understanding global region functionality in Femora