Configure user defined point sources for elastodynamics

configure.ac

@@ -206,6 +206,7 @@ AC_CONFIG_FILES([Makefile
 		libsrc/pylith/topology/Makefile
 		libsrc/pylith/utils/Makefile
 		libsrc/pylith/scales/Makefile
+		libsrc/pylith/sources/Makefile
         libsrc/pylith/testing/Makefile
 		modulesrc/Makefile
 		modulesrc/include/Makefile
@@ -218,6 +219,7 @@ AC_CONFIG_FILES([Makefile
 		modulesrc/mpi/Makefile
 		modulesrc/problems/Makefile
 		modulesrc/scales/Makefile
+		modulesrc/sources/Makefile
 		modulesrc/topology/Makefile
 		modulesrc/utils/Makefile
 		tests/Makefile
@@ -246,6 +248,7 @@ AC_CONFIG_FILES([Makefile
 		tests/mmstests/linearelasticity/nofaults-2d/Makefile
 		tests/mmstests/linearelasticity/nofaults-3d/Makefile
 		tests/mmstests/linearelasticity/faults-2d/Makefile
+		tests/mmstests/linearelasticity/pointsource-2d/Makefile
 		tests/mmstests/incompressibleelasticity/Makefile
 		tests/mmstests/incompressibleelasticity/nofaults-2d/Makefile
 		tests/mmstests/poroelasticity/Makefile
@@ -262,6 +265,7 @@ AC_CONFIG_FILES([Makefile
 		tests/fullscale/linearelasticity/faults-3d/Makefile
 		tests/fullscale/linearelasticity/faults-3d-buried/Makefile
 		tests/fullscale/linearelasticity/greensfns-2d/Makefile
+		tests/fullscale/linearelasticity/pointsource-2d/Makefile
 		tests/fullscale/incompressibleelasticity/Makefile
 		tests/fullscale/incompressibleelasticity/nofaults-2d/Makefile
 		tests/fullscale/incompressibleelasticity/nofaults-3d/Makefile

examples/point-source-2d/Makefile.am

@@ -0,0 +1,20 @@
+# =================================================================================================
+# This code is part of PyLith, developed through the Computational Infrastructure
+# for Geodynamics (https://github.com/geodynamics/pylith).
+#
+# Copyright (c) 2010-2025, University of California, Davis and the PyLith Development Team.
+# All rights reserved.
+#
+# See https://mit-license.org/ and LICENSE.md and for license information. 
+# =================================================================================================
+
+dist_noinst_DATA = \
+	README.md \
+	pylithapp.cfg \
+	step01_explosion.cfg \
+	step02_strikeslip.cfg \
+	step03_multiple.cfg \
+	step04_ramp.cfg \
+	step05_step.cfg
+
+# End of file

examples/point-source-2d/README.md

@@ -0,0 +1,81 @@
+# Point Source Example (2D)
+
+This example demonstrates how to use moment tensor point sources in
+elastodynamic simulations with PyLith. Point sources are useful for
+simulating earthquake sources without explicitly modeling a fault.
+
+## Features
+
+- Point moment tensor sources with user-defined:
+  - Location (x, y coordinates)
+  - Moment tensor components
+  - Magnitude (seismic moment)
+  - Source time function (Ricker wavelet, Gaussian, step, ramp, or custom)
+
+## Simulations
+
+### Step 1: Explosion Source
+
+A simple isotropic (explosion) source demonstrating:
+- Location at the center of the domain
+- Isotropic moment tensor [Mxx=1, Myy=1, Mzz=1, Mxy=0]
+- Ricker wavelet source time function
+
+### Step 2: Strike-slip Source
+
+A double-couple source typical of strike-slip earthquakes:
+- Location offset from center
+- Strike-slip moment tensor [Mxx=1, Myy=-1, Mxy=0]
+- Gaussian source time function
+
+### Step 3: Multiple Sources
+
+Multiple point sources in the same simulation:
+- Two sources at different locations
+- Different source mechanisms
+- Different source time functions
+
+## Running the Simulations
+
+```bash
+# Step 1: Explosion source
+pylith step01_explosion.cfg
+
+# Step 2: Strike-slip source  
+pylith step02_strikeslip.cfg
+
+# Step 3: Multiple sources
+pylith step03_multiple.cfg
+```
+
+## Configuration
+
+The point sources are configured in the `[pylithapp.problem.sources]` section.
+Key parameters include:
+
+- `location`: Source coordinates [x, y] or [x, y, z]
+- `moment_tensor`: Moment tensor components in Voigt notation
+- `magnitude`: Seismic moment (N*m)
+- `time_function`: Source time function type and parameters
+
+See the individual step configuration files for detailed examples.
+
+## Moment Tensor Convention
+
+The moment tensor is specified in Voigt notation:
+- 2D: [Mxx, Myy, Mxy]
+- 3D: [Mxx, Myy, Mzz, Mxy, Mxz, Myz]
+
+Common source mechanisms:
+- Explosion: [1, 1, 1, 0, 0, 0] (isotropic)
+- Strike-slip: [1, -1, 0, 0, 0, 0]
+- Vertical dipole: [0, 0, 1, 0, 0, 0]
+
+## Source Time Functions
+
+Available source time functions:
+- `SourceTimeStep`: Step (Heaviside) function
+- `SourceTimeRamp`: Linear ramp function
+- `SourceTimeGaussian`: Gaussian pulse
+- `SourceTimeRicker`: Ricker wavelet (Mexican hat)
+- `SourceTimeHistory`: User-defined from file

examples/point-source-2d/pylithapp.cfg

@@ -0,0 +1,167 @@
+[pylithapp.metadata]
+# This is not a self-contained simulation configuration file. This
+# file only specifies the general parameters common to the simulations
+# in this directory.
+keywords = [example, 2D, point source, moment tensor, elastodynamics]
+features = [
+    Dynamic simulation,
+    Point source,
+    Moment tensor,
+    Quadrilateral cells,
+    Runge-Kutta time stepping,
+    pylith.meshio.MeshIOCubit,
+    pylith.problems.TimeDependent,
+    pylith.materials.Elasticity,
+    pylith.materials.IsotropicLinearElasticity,
+    pylith.sources.PointMomentTensor,
+    spatialdata.spatialdb.UniformDB,
+    pylith.meshio.DataWriterHDF5
+    ]
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+# Turn on some journals to show progress.
+[pylithapp.journal]
+device = color-console
+
+[pylithapp.journal.info]
+timedependent = 1
+solution = 1
+meshiocubit = 1
+
+[pylithapp.journal.debug]
+timedependent = 0
+
+# ----------------------------------------------------------------------
+# mesh_generator
+# ----------------------------------------------------------------------
+[pylithapp.mesh_generator]
+# Use the same mesh as barwaves-2d example
+# Set the reader to match the type of mesh file.
+reader = pylith.meshio.MeshIOCubit
+reader.filename = ../barwaves-2d/quad_200m.exo
+
+# Set the Cartesian coordinate system.
+reader.coordsys.space_dim = 2
+
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Use the nonlinear solver to verify residual and Jacobian are consistent.
+solver = linear
+formulation = dynamic
+
+# Nondimensionalize problem using wave propagation parameters.
+scales = pylith.scales.NondimElasticDynamic
+scales.mass_density = 2500.0*kg/m**3
+scales.shear_wave_speed = 1.0*km/s
+scales.wave_period = 2.0*s
+
+defaults.quadrature_order = 1
+
+
+# Set the discretization for each of the solution subfields.
+solution = pylith.problems.SolnDispVel
+
+solution.subfields.displacement.basis_order = 1
+solution.subfields.velocity.basis_order = 1
+
+solution_observers = [domain]
+
+start_time = -0.2*s
+end_time = 10.0*s
+initial_dt = 0.1*s
+
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Create an array of one material
+materials = [elastic]
+
+# We use the default material (elasticity) and rheology
+# (isotropic, linearly elastic).
+
+[pylithapp.problem.materials.elastic]
+# label_value must match the block values in the Cubit Exodus file.
+label_value = 1
+
+# We will use uniform material properties, so we use the UniformDB
+# spatial database.
+db_auxiliary_field = spatialdata.spatialdb.UniformDB
+db_auxiliary_field.description = Elastic properties
+db_auxiliary_field.values = [density, vs, vp]
+db_auxiliary_field.data = [2500*kg/m**3, 1.0*km/s, 1.732*km/s]
+
+# Set the discretization of the material auxiliary fields (properties).
+# We have uniform material properties, so we can use a basis order of 0.
+auxiliary_subfields.density.basis_order = 0
+bulk_rheology.auxiliary_subfields.bulk_modulus.basis_order = 0
+bulk_rheology.auxiliary_subfields.shear_modulus.basis_order = 0
+
+
+# ----------------------------------------------------------------------
+# boundary conditions
+# ----------------------------------------------------------------------
+# Set boundary conditions to absorb waves at domain boundaries
+[pylithapp.problem]
+bc = [bc_xneg, bc_xpos, bc_yneg, bc_ypos]
+bc.bc_xneg = pylith.bc.AbsorbingDampers
+bc.bc_xpos = pylith.bc.AbsorbingDampers
+bc.bc_yneg = pylith.bc.AbsorbingDampers
+bc.bc_ypos = pylith.bc.AbsorbingDampers
+
+[pylithapp.problem.bc.bc_xneg]
+label = boundary_xneg
+db_auxiliary_field = spatialdata.spatialdb.UniformDB
+db_auxiliary_field.description = Absorbing BC -x
+db_auxiliary_field.values = [density, vs, vp]
+db_auxiliary_field.data = [2500*kg/m**3, 1.0*km/s, 1.732*km/s]
+
+[pylithapp.problem.bc.bc_xpos]
+label = boundary_xpos
+db_auxiliary_field = spatialdata.spatialdb.UniformDB
+db_auxiliary_field.description = Absorbing BC +x
+db_auxiliary_field.values = [density, vs, vp]
+db_auxiliary_field.data = [2500*kg/m**3, 1.0*km/s, 1.732*km/s]
+
+[pylithapp.problem.bc.bc_yneg]
+label = boundary_yneg
+db_auxiliary_field = spatialdata.spatialdb.UniformDB
+db_auxiliary_field.description = Absorbing BC -y
+db_auxiliary_field.values = [density, vs, vp]
+db_auxiliary_field.data = [2500*kg/m**3, 1.0*km/s, 1.732*km/s]
+
+[pylithapp.problem.bc.bc_ypos]
+label = boundary_ypos
+db_auxiliary_field = spatialdata.spatialdb.UniformDB
+db_auxiliary_field.description = Absorbing BC +y
+db_auxiliary_field.values = [density, vs, vp]
+db_auxiliary_field.data = [2500*kg/m**3, 1.0*km/s, 1.732*km/s]
+
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[pylithapp.petsc]
+ts_type = rk
+ts_rk_type = 3bs
+ts_adapt_dt_max = 0.05
+
+ksp_rtol = 1.0e-8
+ksp_atol = 1.0e-12
+ksp_max_it = 30
+ksp_gmres_restart = 50
+ksp_error_if_not_converged = true
+
+snes_rtol = 1.0e-10
+snes_atol = 1.0e-10
+snes_error_if_not_converged = true
+
+# Monitors for debugging
+ts_monitor = true
+
+# End of file

examples/point-source-2d/step01_explosion.cfg

@@ -0,0 +1,68 @@
+[pylithapp.metadata]
+#  y
+#  ^
+#  |
+#   --> x
+#
+# Domain: -4 km <= x <= 4 km, -4 km <= y <= 4 km
+#
+# Isotropic (explosion) point source at center of domain.
+# Absorbing boundary conditions on all sides.
+#
+base = [pylithapp.cfg]
+description = Explosion point source with Ricker wavelet time function.
+authors = [PyLith Development Team]
+keywords = [explosion, isotropic source, Ricker wavelet, point source]
+arguments = [step01_explosion.cfg]
+version = 1.0.0
+pylith_version = [>=4.0, <6.0]
+
+features = [
+    pylith.sources.PointMomentTensor,
+    pylith.sources.SourceTimeRicker
+    ]
+
+
+# ----------------------------------------------------------------------
+# output
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Set the name of the problem that will be used to construct the
+# output filenames. The default directory for output is 'output'.
+defaults.name = step01_explosion
+
+
+# ----------------------------------------------------------------------
+# point sources
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Define point sources for the simulation
+sources = [explosion]
+
+# Configure the explosion source
+[pylithapp.problem.sources.explosion]
+# Location at center of domain
+location = [0.0*km, 0.0*km]
+
+# Isotropic (explosion) moment tensor
+# 2D format: [Mxx, Myy, Mxy]
+# Explosion has equal diagonal components and zero shear
+moment_tensor = [1.0, 1.0, 0.0]
+
+# Seismic moment magnitude (in N*m)
+# M0 = 1e15 N*m corresponds to approximately Mw 4.0
+magnitude = 1.0e+15*N*m
+
+# Source time function: Ricker wavelet
+time_function = pylith.sources.SourceTimeRicker
+time_function.origin_time = 0.0*s
+time_function.peak_frequency = 2.0*Hz
+time_function.delay = 1.0*s
+time_function.magnitude = 1.0
+
+# Characteristic length for source regularization
+# Should be similar to element size
+characteristic_length = 200.0*m
+
+
+# End of file