Merge pull request #667 from baagaard-usgs/feature-output-traction-ch…

…ange

Add output of fault traction changes and boundary and fault orientation information

CHANGES.md

@@ -1,16 +1,27 @@
 See <https://github.com/geodynamics/pylith/commits/main> for the complete log of changes made to PyLith.
 
+:::{note}
+Starting with v3.0.0, we strictly follow the [semantic versioning guidelines](https://semver.org/).
+The version numbers are in the form `MAJOR.MINOR.PATCH`, where major releases indicate changes to the public API (parameters), minor releases indicate new functionality that is backward compatible, and patch releases indicate backward compatible bug fixes.
+:::
+
 ## Version 4.0.0
 
+### Changes in user parameters
+
+* Changed name of fault Lagrange multiplier field for solution component in Python from `lagrange_fault` to `lagrange_multiplier_fault` to match name of solution field in C++.
 * Removed support for importing meshes from LaGriT.
-* Add output of change in fault tractions for prescribed slip.
+
+### Other changes
+
+* Change in fault tractions are now included in the fault `data_fields` for prescribed slip.
+* Fault and boundary orientation directions are now included in the `info_fields` for simulation output.
 * State variables are now included in the default `data_fields` for simulation output.
 * The default solver settings use the PETSc proper orthogonal decomposition (POD) methodology for initial guess of solutions to improve convergence.
-* Changed name of fault Lagrange multiplier field for solution component in Python from `lagrange_fault` to `lagrange_multiplier_fault` to match name of solution field in C++.
 * Add demonstration of `pylith_powerlaw_gendb` in Step 8 of `examples/reverse-2d`.
 * Add demonstration of using poroelasticity with porosity as a state variable to `examples/magma-2d`.
 * Switched from CppUnit to Catch2 as the C++ testing framework.
-* Update to PETSc 3.19.5
+* Update to PETSc 3.20.2
 * Improve integration with VSCode for testing and debugging (see Developer Guide)
 * Bug fixes
   * Fix errors in KinSrcTimeHistory.py

configure.ac

@@ -109,7 +109,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])])
 
 dnl PETSC
 AC_LANG(C)
-CIT_PATH_PETSC([3.19.5])
+CIT_PATH_PETSC([3.20.2])
 CIT_HEADER_PETSC
 CIT_CHECK_LIB_PETSC
 

docs/intro/development-plan.md

@@ -6,48 +6,53 @@ Future implementation of features is guided by several target applications, incl
 * Inversion of geodetic data for slow slip events, fault creep, and long-term fault slip rates.
 * Quasistatic and dynamic modeling of fluids and faulting.
 
-## Version 3.0.4 (September 2023)
+:::{note}
+Because we strictly follow the [semantic versioning guidelines](https://semver.org/), a minor release may get promoted to a major releases if we make changes to the public API (parameters).
+This can happen if realize that we should modify the parameters to improve maintainability or prepare for future changes.
+:::
 
-Updates for examples and documentation along with bugfixes.
+## Version 4.0.0 (December 2023)
 
 * Analytical function spatial database ![easy](images/easy.png)[100%]
 * Use `pod` initial guess to improve convergence ![easy](images/easy.png) [100%]
-* Dynamic prescribed slip with diagonal Jacobian for explicit part of IMEX formulation ![expert](images/expert.png)[75%]
-* Output of fault tractions ![expert](images/expert.png) [75%]
-* Update coordinates with solution ![intermediate](images/intermediate.png) [50%]
-* Convert from CppUnit to Catch2 ![easy](images/easy.png) [40%]
+* Output of fault tractions ![expert](images/expert.png) [100%]
+* Output of boundary and fault orientation ![export](images/expert.png) [100%]
+* Convert from CppUnit to Catch2 ![easy](images/easy.png) [100%]
+
+## Version 4.1.0 (March 2024)
+
+* Output of fault rupture auxiliary subfields ![intermediate](images/intermediate.png) [0%]
+* Dynamic prescribed slip with diagonal Jacobian for explicit part of IMEX formulation ![expert](images/expert.png) [75%]
+* Parallel mesh loading ![expert](images/expert.png) [15%]
+* Add 2D and 3D examples for crustal faults with complex fault geometry ![easy](images/easy.png) [50%]
 * Better preconditioners ![expert](images/expert.png) [25%]
   * elasticity with fault
   * incompressible elasticity
   * poroelasticity
-* Finish updating `examples/subduction-3d` ![intermediate](images/intermediate.png) [20%]
-* Add `examples/barwaves-2d` ![expert](images/expert.png) [15%]
-* Parallel mesh loading ![expert](images/expert.png) [15%]
-* Elasticity with self-gravitation ![intermediate](images/intermediate.png) [0%]
-* Add 2D and 3D examples for crustal faults with complex fault geometry ![easy](images/easy.png) [0%]
+* Update VTK output to use `vtu` files rather than legacy `vtk` files ![easy](images/easy.png) [0%]
+* Finish updating `examples/subduction-3d` ![intermediate](images/intermediate.png) [40%]
 
-## Version 3.1.0 (December 2023)
+## Version 5.0.0 (June 2024)
 
-* Spontaneous rupture for quasistatic simulations ![expert](images/expert.png) [20%]
+* Spontaneous rupture for quasistatic and dynamic simulations ![expert](images/expert.png) [20%]
+* Convert from SWIG to pybind11 ![intermediate](images/intermediate.png) [0%]
+* Reimplementation of Drucker-Prager elastoplastic bulk rheology ![intermediate](images/intermediate.png) [0%]
 * Add support for GeoModelGrids implementation of spatial databases for 3D seismic velocity models. ![intermediate](images/intermediate.png) [0%]
 * Improve robustness of HDF5 output by opening/closing at each time step ![easy](images/easy.png)[0%]
-* Dirichlet boundary conditions with constraints on normal and tangential components. ![difficult](images/difficult.png) [0%]
-* Additional minor cleanup of code internals to improve maintainability.
-  * Refactor auxiliary field and derived field output (initial/timestep/final)
-  * Output of fault rupture auxiliary fields
-  * VTK output (vtk -> vtu)
-* Reimplementation of Drucker-Prager elastoplastic bulk rheology ![intermediate](images/intermediate.png) [0%]
-* Convert from SWIG to pybind11 ![intermediate](images/intermediate.png) [0%]
+* Add `examples/barwaves-2d` ![expert](images/expert.png) [15%]
+* Update coordinates with solution ![intermediate](images/intermediate.png) [50%]
 
-## Version 4.0 (June 2024)
+## Version 6.0.0 (TBD)
 
-* Spontaneous rupture for dynamic simulations ![expert](images/expert.png) [10%]
 * Migrate examples to Jupyter notebooks ![intermediate](images/intermediate.png)
 * Update to current version of Pyre ![difficult](images/difficult.png)
 * More flexible specification of time-dependent boundary conditions. ![difficult](images/difficult.png) [0%]
+* Dirichlet boundary conditions with constraints on normal and tangential components. ![difficult](images/difficult.png) [0%]
 * Integration with libCEED for fast high order residual evaluation ![expert](images/expert.png)\
   Contribution led by Jed Brown.
 * Add ability to output residual field during nonlinear solve for debugging ![easy](images/easy.png) [0%]
+* Elasticity with self-gravitation ![intermediate](images/intermediate.png) [0%]
+
 
 ## Features for Future Releases
 
@@ -59,7 +64,7 @@ Updates for examples and documentation along with bugfixes.
 * Adaptive mesh refinement ![expert](images/expert.png)
 * Line/point fluid sources in poroelasticity ![expert](images/expert.png) [20%]
 * Consolidate HDF5 output into a single file ![difficult](images/difficult.png)
-* Drucker-Prager bulk rheology with relaxation to yield surface ![intermediate](images/intermediate.png) 
+* Drucker-Prager bulk rheology with relaxation to yield surface ![intermediate](images/intermediate.png)
 * Drucker-Prager bulk rheology with strain hardening/softening  ![intermediate](images/intermediate.png)
 * Adjoint for data assimilation ![difficult](images/difficult.png)
 * Fault with both prescribed slip and spontaneous rupture ![difficult](images/difficult.png)\

docs/user/physics/bc/index.md

@@ -37,6 +37,11 @@ bc.y_pos = pylith.bc.AbsorbingDampers
 bc.z_neg = pylith.bc.NeumannTimeDependent
 ```
 
+## Diagnostic Information
+
+The diagnostic information includes the outward normal direction (`normal_dir`) and the two tangential directions (`horizontal_tangential_dir` and `vertical_tangential_dir`).
+The default basis order for discretizing these directions is 1, so these produce `vertex_fields` as opposed to `cell_fields` (basis order of 0).
+
 ## Boundary Condition Implementations
 
 :::{toctree}

docs/user/physics/faults/index.md

@@ -68,7 +68,12 @@ The default value of [0, 0, 1] is appropriate for most 3D problems.
 
 By default the output observers write both diagnostic information (for example, fault orientation directions) and the slip at each time step.
 The fault coordinate system is shown in {numref}`fig:fault:slip:motions`.
-The vectors in the fault coordinate system can be transformed to the global coordinate system using the direction vectors in the diagnostic output.
+The vectors in the fault coordinate system can be transformed to the global coordinate system using the direction vectors in the diagnostic output ({ref}`sec-user-output-observers`).
+
+:::{important}
+The normal direction is chosen based on how the cells are split to create cohesive cells.
+If the normal direction contains a positive z component, then the directions conform to traditional seismologic conventions (along strike and up dip); however, if the normal direction contains a negative z component, then the directions correspond to along strike and down dip directions.
+:::
 
 :::{toctree}
 prescribed-slip.md

docs/user/physics/faults/prescribed-slip.md

@@ -284,3 +284,18 @@ The slip and slip initiation times are specified using spatial databases, so the
 :::{seealso}
 See [`KinSrcTimeHistory` Component](../components/faults/../../../components/faults/KinSrcTimeHistory.md) for the Pyre properties and facilities and configuration examples.
 :::
+
+## Output
+
+The derived subfield available for prescribed slip is the change in tractions on the fault surface.
+
+:::{tip}
+To compute the change in tractions on a locked fault, prescribe zero slip on the fault.
+:::
+
+```{table} Derived subfields that are available for output for prescribed slip.
+:name: tab:fault:prescribed_slip:derived:subfields
+|      Subfield |   Components            |
+|:--------------|:-----------------------|
+| `traction_change` | normal, along strike, up dip |
+```

docs/user/physics/faults/slip-impulses.md

@@ -25,3 +25,14 @@ Impulses will be applied at any point on the fault with a slip component greater
 :class: seealso
 See [`FaultCohesiveImpulses` Component](../../components/faults/FaultCohesiveImpulses.md)  for the Pyre properties and facilities and configuration examples.
 :::
+
+## Output
+
+The derived subfield available for Green's functions is the change in tractions on the fault surface.
+
+```{table} Derived subfields that are available for output for Green's functions.
+:name: tab:fault:prescribed_slip:derived:subfields
+|      Subfield |   Components            |
+|:--------------|:-----------------------|
+| `traction_change` | normal, along strike, up dip |
+```

docs/user/problems/output.md

@@ -36,8 +36,9 @@ The monitor calculates the percent completed based on the number of steps (e.g.,
 ## Observers
 
 Observer objects manage transferring the solution and state variables to other objects, including output and external software.
-Currently, PyLith only includes observers for output.
+Currently, the only type of observers implemented in PyLith are ones that produce output.
 
+(sec-user-output-observers)=
 ### Output Observers
 
 PyLith currently supports output to HDF5/Xdmf and VTK files, which can be imported directly into a number of visualization tools, such as ParaView, Visit, and MayaVi.
@@ -55,6 +56,27 @@ Fields with a basis order of 0 are kept at a basis order of 0 when output.
 Fields with a basis order of 1 or more can be output with a basis order of 0 or 1.
 :::
 
+The output observers produce files with either diagnostic information (`info` files) or solution and state variable information (`data` files).
+The default behavior is the files include all available information.
+The `info` files include the auxiliary subfields at the beginning of the simulation along with  surface orientation information for faults and boundary conditions.
+The `data` files include all solution subfields, state variables, and derived fields (fields computed from the solution, such as Cauchy stress and strain).
+
+For boundary conditions the orientation information is provided in terms of x, y, and z components of the unit vectors for the surface normal and tangential directions.
+In 3D the "vertical" tangential direction is the cross product of the surface normal and horizontal tangential direction; it is in the +z direction for a vertical boundary.
+In the case of the horizontal boundary, the horizontal tangential direction is in the +x direction and the "vertical" tangential direction is in the +y direction.
+For a fault surface the horizontal tangential direction generally corresponds to the along-strike direction and the "vertical" tangential direction generally corresponds to the up-dip direction; the exception is a 2D simulation for a vertical cross-section in which the "horizontal" tangential direction corresponds to the dip direction.
+
+The orientation information is useful for transforming from components written in terms of a surface (for example, left lateral, reverse, and opening fault tractions), into the model coordinate system (for example, the global Cartesian coordinate system).
+Given unit vector components $n_x$, $n_y$, $n_z$ (normal direction), $h_x$, $h_y$, $h_z$ (horizontal tangential direction or along-strike direction), and $v_x$, $v_y$, $v_z$ (vertical tangential direction or up-dip direction), we can transform a vector in the boundary (or fault) coordinate system ($a_n$, $a_h$, $a_v$) into the global coordinate system using
+
+\begin{align}
+a_x &= a_n n_x + a_h h_x + a_v v_x \\
+a_y &= a_n n_y + a_h h_y + a_v v_y \\
+a_z &= a_n n_z + a_h h_z + a_v v_z
+\end{align}
+
+This transformation is useful for plotting fault tractions and slip in 3D visualization tools such as ParaView that require vectors in the model coordinate system.
+
 (sec-user-solution-observers)=
 ### Solution Observers
 

libsrc/pylith/Makefile.am

@@ -33,6 +33,7 @@ lib_LTLIBRARIES = libpylith.la
 
 libpylith_la_SOURCES = \
 	bc/BoundaryCondition.cc \
+	bc/DiagnosticFieldFactory.cc \
 	bc/TimeDependentAuxiliaryFactory.cc \
 	bc/DirichletTimeDependent.cc \
 	bc/DirichletUserFn.cc \
@@ -43,7 +44,9 @@ libpylith_la_SOURCES = \
 	faults/FaultCohesive.cc \
 	faults/FaultCohesiveKin.cc \
 	faults/FaultCohesiveImpulses.cc \
-	faults/AuxiliaryFactoryKinematic.cc \
+	faults/AuxiliaryFieldFactory.cc \
+	faults/DiagnosticFieldFactory.cc \
+	faults/DerivedFieldFactory.cc \
 	faults/KinSrc.cc \
 	faults/KinSrcStep.cc \
 	faults/KinSrcRamp.cc \
@@ -120,6 +123,7 @@ libpylith_la_SOURCES = \
 	problems/ObserverSoln.cc \
 	problems/ObserversSoln.cc \
 	problems/Physics.cc \
+	problems/Observer.cc \
 	problems/ObserverPhysics.cc \
 	problems/ObserversPhysics.cc \
 	problems/InitialCondition.cc \

libsrc/pylith/bc/AbsorbingDampers.cc

@@ -131,6 +131,7 @@ pylith::bc::AbsorbingDampers::createIntegrator(const pylith::topology::Field& so
     integrator->setLabelValue(getLabelValue());
 
     _AbsorbingDampers::setKernelsResidual(integrator, *this, solution);
+    BoundaryCondition::_setKernelsDiagnosticField(integrator, solution);
 
     PYLITH_METHOD_RETURN(integrator);
 } // createIntegrator
@@ -188,18 +189,6 @@ pylith::bc::AbsorbingDampers::createAuxiliaryField(const pylith::topology::Field
 } // createAuxiliaryField
 
 
-// ---------------------------------------------------------------------------------------------------------------------
-// Create derived field.
-pylith::topology::Field*
-pylith::bc::AbsorbingDampers::createDerivedField(const pylith::topology::Field& solution,
-                                                 const pylith::topology::Mesh& domainMesh) {
-    PYLITH_METHOD_BEGIN;
-    PYLITH_COMPONENT_DEBUG("createDerivedField(solution="<<solution.getLabel()<<", domainMesh=)"<<typeid(domainMesh).name()<<") empty method");
-
-    PYLITH_METHOD_RETURN(NULL);
-} // createDerivedField
-
-
 // ---------------------------------------------------------------------------------------------------------------------
 // Get auxiliary field factory associated with physics.
 pylith::feassemble::AuxiliaryFactory*
@@ -208,25 +197,6 @@ pylith::bc::AbsorbingDampers::_getAuxiliaryFactory(void) {
 } // _getAuxiliaryFactory
 
 
-// ---------------------------------------------------------------------------------------------------------------------
-// Update kernel constants.
-void
-pylith::bc::AbsorbingDampers::_updateKernelConstants(const PylithReal dt) {
-    PYLITH_METHOD_BEGIN;
-    PYLITH_COMPONENT_DEBUG("_setKernelConstants(dt="<<dt<<")");
-
-    if (6 != _kernelConstants.size()) { _kernelConstants.resize(6);}
-    _kernelConstants[0] = _refDir1[0];
-    _kernelConstants[1] = _refDir1[1];
-    _kernelConstants[2] = _refDir1[2];
-    _kernelConstants[3] = _refDir2[0];
-    _kernelConstants[4] = _refDir2[1];
-    _kernelConstants[5] = _refDir2[2];
-
-    PYLITH_METHOD_END;
-} // _updateKernelConstants
-
-
 // ---------------------------------------------------------------------------------------------------------------------
 // Set kernels for residual.
 void

libsrc/pylith/bc/AbsorbingDampers.hh

@@ -73,16 +73,6 @@ public:
     pylith::topology::Field* createAuxiliaryField(const pylith::topology::Field& solution,
                                                   const pylith::topology::Mesh& domainMesh);
 
-    /** Create derived field.
-     *
-     * @param[in] solution Solution field.
-     * @param[in\ domainMesh Finite-element mesh associated with integration domain.
-     *
-     * @returns Derived field if applicable, otherwise NULL.
-     */
-    pylith::topology::Field* createDerivedField(const pylith::topology::Field& solution,
-                                                const pylith::topology::Mesh& domainMesh);
-
     // PROTECTED METHODS ///////////////////////////////////////////////////////////////////////////////////////////////
 protected:
 
@@ -92,18 +82,10 @@ protected:
      */
     pylith::feassemble::AuxiliaryFactory* _getAuxiliaryFactory(void);
 
-    /** Update kernel constants.
-     *
-     * @param[in] dt Current time step.
-     */
-    void _updateKernelConstants(const PylithReal dt);
-
     // PRIVATE MEMBERS /////////////////////////////////////////////////////////////////////////////////////////////////
 private:
 
     pylith::bc::AbsorbingDampersAuxiliaryFactory* _auxiliaryFactory; ///< Factory for auxiliary subfields.
-    PylithReal _refDir1[3]; ///< First choice reference direction used to compute boundary tangential directions.
-    PylithReal _refDir2[3]; ///< Second choice reference direction used to compute boundary tangential directions.
 
     // NOT IMPLEMENTED /////////////////////////////////////////////////////////////////////////////////////////////////
 private: