Point source merge

.vscode/settings.json

@@ -10,7 +10,9 @@
     },
     "files.associations": {
         "*.icc": "cpp",
-        "iosfwd": "cpp"
+        "iosfwd": "cpp",
+        "*.tcc": "cpp",
+        "fstream": "cpp"
     },
     "python.formatting.autopep8Args": [
         "-max-line-length",

configure.ac

@@ -199,6 +199,7 @@ AC_CONFIG_FILES([Makefile
 		libsrc/pylith/feassemble/Makefile
 		libsrc/pylith/fekernels/Makefile
 		libsrc/pylith/faults/Makefile
+		libsrc/pylith/sources/Makefile
 		libsrc/pylith/friction/Makefile
 		libsrc/pylith/materials/Makefile
 		libsrc/pylith/meshio/Makefile
@@ -210,6 +211,7 @@ AC_CONFIG_FILES([Makefile
 		modulesrc/include/Makefile
 		modulesrc/bc/Makefile
 		modulesrc/faults/Makefile
+		modulesrc/sources/Makefile
 		modulesrc/feassemble/Makefile
 		modulesrc/friction/Makefile
 		modulesrc/materials/Makefile

docs/user/governingeqns/elasticity-infstrain/dynamic.md

@@ -9,7 +9,7 @@ We introduce the velocity as a unknown, $\vec{v}=\frac{\partial u}{\partial t}$,
 % Displacement-velocity
 \frac{\partial \vec{u}}{\partial t} &= \vec{v} \text{ in } \Omega, \\
 % Elasticity
-\rho(\vec{x}) \frac{\partial\vec{v}}{\partial t} &= \vec{f}(\vec{x},t) + \boldsymbol{\nabla} \cdot \boldsymbol{\sigma}(\vec{u}) \text{ in } \Omega, \\
+\rho(\vec{x}) \frac{\partial\vec{v}}{\partial t} &= \vec{f}(\vec{x},t) + \boldsymbol{\nabla} \cdot \boldsymbol{\sigma}(\vec{u}) + \boldsymbol{\nabla} \cdot {\bf M} \text{ in } \Omega, \\
 % Neumann
 \boldsymbol{\sigma} \cdot \vec{n} &= \vec{\tau}(\vec{x},t) \text{ on } \Gamma_\tau, \\
 % Dirichlet

examples/point-source/impulse_gen.py

@@ -0,0 +1,29 @@
+import numpy as np
+
+def generate_ramp_stf(duration=0.5, dt=0.01, max_amplitude=1.0):
+    time = np.arange(0, duration + dt, dt)
+    amplitude = max_amplitude - (max_amplitude / duration) * time
+    return time, amplitude
+
+# Parameters
+duration = 0.10
+dt = 0.001
+max_amp = 1.0
+time, amplitude = generate_ramp_stf(duration, dt, max_amp)
+
+# Write in the same format as impulse.timedb
+with open("ramp.timedb", "w") as f:
+    f.write("// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting)\n")
+    f.write("//\n")
+    f.write("// Ramp source time function with linear increase to max amplitude.\n")
+    f.write("#TIME HISTORY ascii\n")
+    f.write("TimeHistory {\n")
+    f.write(f"  num-points = {len(time)+1} // number of points in time history\n")
+    f.write("  time-units = second // units for time\n")
+    f.write("}\n")
+    for t, a in zip(time, amplitude):
+        f.write(f"{t:10.4f}   {a:10.4f}\n")
+    t_final = 999.000
+    a_final = 0.000
+    f.write(f"{t_final:10.4f}   {a_final:10.4f}\n")
+

examples/point-source/mesh_quad.jou

@@ -0,0 +1,30 @@
+# ----------------------------------------------------------------------
+# Generate geometry
+# ----------------------------------------------------------------------
+playback 'geometry.jou'
+
+# ----------------------------------------------------------------------
+# Set discretization size
+# ----------------------------------------------------------------------
+${dx=120.0*m}
+surface all size {dx}
+
+
+# ----------------------------------------------------------------------
+# Generate the mesh
+# ----------------------------------------------------------------------
+surface all scheme submap
+mesh surface all
+
+# ----------------------------------------------------------------------
+# Mark entities for boundary conditions, etc.
+# ----------------------------------------------------------------------
+playback 'bc.jou'
+
+# ----------------------------------------------------------------------
+# Export exodus file
+# ----------------------------------------------------------------------
+export mesh "mesh_quad.exo" dimension 2 overwrite
+
+
+

examples/point-source/pylithapp.cfg

@@ -0,0 +1,124 @@
+[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]
+features = [
+    Dynamic simulation,
+    Implicit-explicit formulation,
+    Quadrilateral cells,
+    Runge-Kutta time stepping,
+    pylith.meshio.MeshIOCubit,
+    pylith.problems.TimeDependent,
+    pylith.materials.Elasticity,
+    pylith.materials.IsotropicLinearElasticity,
+    spatialdata.spatialdb.UniformDB,
+    pylith.meshio.DataWriterHDF5,
+    spatialdata.units.NondimElasticDynamic
+    ]
+
+# ----------------------------------------------------------------------
+# journal
+# ----------------------------------------------------------------------
+# Turn on some journals to show progress.
+[pylithapp.journal]
+#device = color-console
+
+[pylithapp.journal.info]
+timedependent = 1
+solution = 1
+meshiocubit = 1
+pylithapp = 1
+
+[pylithapp.journal.debug]
+timedependent = 0
+
+
+
+# ----------------------------------------------------------------------
+# problem
+# ----------------------------------------------------------------------
+[pylithapp.problem]
+# Use the nonlinear solver to verify residual and Jacobian are consistent.
+solver = linear
+formulation = dynamic
+
+# # Nondimensionalize problem using wave propagation parameters.
+# normalizer = spatialdata.units.NondimElasticDynamic
+# normalizer.mass_density = 2500.0*kg/m**3
+# normalizer.shear_wave_speed = 1.0*km/s
+# normalizer.wave_period = 0.5*s
+
+defaults.quadrature_order = 1
+
+
+# Set the discretization for each of the solution subfields.
+#
+# For a quastistatic simulation with a fault, we have two solution fields:
+# (1) displacement and (2) Lagrange multiplier. We use a predefined containter
+# to create a solution field with these two subfields.
+solution = pylith.problems.SolnDispVel
+
+solution.subfields.displacement.basis_order = 1
+solution.subfields.velocity.basis_order = 1
+
+#solution_observers = [domain]
+#solution_observers = []
+#solution_observers.domain.trigger.num_skip = 10
+
+
+start_time = 0.0*s
+end_time = 0.025*s
+initial_dt = 0.001*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]
+# id must match the values in the mesh material-ids.
+description = Elastic material
+label_value = 1
+observers = []
+#observers.observer.trigger.num_skip = 10
+# 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 = [1800*kg/m**3, 1.3*km/s, 3.0*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
+
+# ----------------------------------------------------------------------
+# PETSc
+# ----------------------------------------------------------------------
+[pylithapp.petsc]
+ts_type = rk
+ts_rk_type = 3bs
+ts_adapt_dt_max = 0.1
+
+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/ramp.timedb

@@ -0,0 +1,110 @@
+// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting)
+//
+// Ramp source time function with linear increase to max amplitude.
+#TIME HISTORY ascii
+TimeHistory {
+  num-points = 102 // number of points in time history
+  time-units = second // units for time
+}
+    0.0000       1.0000
+    0.0010       0.9900
+    0.0020       0.9800
+    0.0030       0.9700
+    0.0040       0.9600
+    0.0050       0.9500
+    0.0060       0.9400
+    0.0070       0.9300
+    0.0080       0.9200
+    0.0090       0.9100
+    0.0100       0.9000
+    0.0110       0.8900
+    0.0120       0.8800
+    0.0130       0.8700
+    0.0140       0.8600
+    0.0150       0.8500
+    0.0160       0.8400
+    0.0170       0.8300
+    0.0180       0.8200
+    0.0190       0.8100
+    0.0200       0.8000
+    0.0210       0.7900
+    0.0220       0.7800
+    0.0230       0.7700
+    0.0240       0.7600
+    0.0250       0.7500
+    0.0260       0.7400
+    0.0270       0.7300
+    0.0280       0.7200
+    0.0290       0.7100
+    0.0300       0.7000
+    0.0310       0.6900
+    0.0320       0.6800
+    0.0330       0.6700
+    0.0340       0.6600
+    0.0350       0.6500
+    0.0360       0.6400
+    0.0370       0.6300
+    0.0380       0.6200
+    0.0390       0.6100
+    0.0400       0.6000
+    0.0410       0.5900
+    0.0420       0.5800
+    0.0430       0.5700
+    0.0440       0.5600
+    0.0450       0.5500
+    0.0460       0.5400
+    0.0470       0.5300
+    0.0480       0.5200
+    0.0490       0.5100
+    0.0500       0.5000
+    0.0510       0.4900
+    0.0520       0.4800
+    0.0530       0.4700
+    0.0540       0.4600
+    0.0550       0.4500
+    0.0560       0.4400
+    0.0570       0.4300
+    0.0580       0.4200
+    0.0590       0.4100
+    0.0600       0.4000
+    0.0610       0.3900
+    0.0620       0.3800
+    0.0630       0.3700
+    0.0640       0.3600
+    0.0650       0.3500
+    0.0660       0.3400
+    0.0670       0.3300
+    0.0680       0.3200
+    0.0690       0.3100
+    0.0700       0.3000
+    0.0710       0.2900
+    0.0720       0.2800
+    0.0730       0.2700
+    0.0740       0.2600
+    0.0750       0.2500
+    0.0760       0.2400
+    0.0770       0.2300
+    0.0780       0.2200
+    0.0790       0.2100
+    0.0800       0.2000
+    0.0810       0.1900
+    0.0820       0.1800
+    0.0830       0.1700
+    0.0840       0.1600
+    0.0850       0.1500
+    0.0860       0.1400
+    0.0870       0.1300
+    0.0880       0.1200
+    0.0890       0.1100
+    0.0900       0.1000
+    0.0910       0.0900
+    0.0920       0.0800
+    0.0930       0.0700
+    0.0940       0.0600
+    0.0950       0.0500
+    0.0960       0.0400
+    0.0970       0.0300
+    0.0980       0.0200
+    0.0990       0.0100
+    0.1000       0.0000
+  999.0000       0.0000

examples/point-source/source_sites.txt

@@ -0,0 +1,2 @@
+#Site	X	    Y   
+s000	5000.00	5000.00