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

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:

\[[C] = \alpha_M [M] + \beta_K [K]\]

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:

\[ \begin{align}\begin{aligned}\omega_1 = 2\pi f_1, \quad \omega_2 = 2\pi f_2\\\alpha_M = \frac{2\omega_1\omega_2}{\omega_1 + \omega_2}\xi\\\beta_K = \frac{2}{\omega_1 + \omega_2}\xi\end{aligned}\end{align} \]

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:

\[\xi(\omega) = \frac{\alpha_M}{2\omega} + \frac{\beta_K\omega}{2}\]

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 part
- Without 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 waves
thickness: 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:

Example 2 - Multi-layer Soil Model with Absorbing Boundaries

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