Merge branch 'geodynamics:main' into point_source

CHANGES.md

@@ -3,8 +3,15 @@ See <https://github.com/geodynamics/pylith/commits/main> for the complete log of
 ## Version 3.1.0
 
 * Removed support for importing meshes from LaGriT.
+* Add output of change in fault tractions for prescribed slip.
+* By default use PETSc proper orthogonal decomposition (POD) methodology for initial guess of solutions to improve convergence.
+* Add demonstration of `pylith_powerlaw_gendb` in Step 8 of `examples/reverse-2d`.
 * Switched from CppUnit to Catch2 as the C++ testing framework.
+* Update to PETSc 3.19.5
+* Improve integration with VSCode for testing and debugging (see Developer Guide)
 * Bug fixes
+  * Fix errors in KinSrcTimeHistory.py
+  * Fix creation of PETSc label for edges when importing Gmsh files. This fixes creation of faults with buried edges for 3D meshes imported from Gmsh.
 
 ## Version 3.0.3
 

applications/pylith_powerlaw_gendb

@@ -28,7 +28,8 @@ from pythia.pyre.applications.Script import Script as Application
 class PowerLawApp(Application):
     """Pyre application to generate a spatial database with power-law parameters for PyLith.
     """
-
+    TENSOR_COMPONENTS = ["xx", "yy", "zz", "xy", "yz", "xz"]
+
     import pythia.pyre.inventory
     from pythia.pyre.units.pressure import MPa
     from pythia.pyre.units.time import s
@@ -40,11 +41,11 @@ class PowerLawApp(Application):
     refSelection.meta['tip'] = "Indicates whether reference stress or " \
         "reference strain rate is provided as input."
 
-    refStress = pythia.pyre.inventory.dimensional("reference_stress", default=1.0 * MPa)
-    refStress.meta['tip'] = "Reference stress value."
+    powerlawRefStress = pythia.pyre.inventory.dimensional("reference_stress", default=1.0 * MPa)
+    powerlawRefStress.meta['tip'] = "Reference stress value."
 
-    refStrainRate = pythia.pyre.inventory.dimensional("reference_strain_rate", default=1.0e-6 / s)
-    refStrainRate.meta['tip'] = "Reference strain rate value."
+    pwerlawRefStrainRate = pythia.pyre.inventory.dimensional("reference_strain_rate", default=1.0e-6 / s)
+    pwerlawRefStrainRate.meta['tip'] = "Reference strain rate value."
 
     from spatialdata.spatialdb.SimpleDB import SimpleDB
 
@@ -60,6 +61,9 @@ class PowerLawApp(Application):
     dbAe = pythia.pyre.inventory.facility("db_powerlaw_coefficient", family="spatial_database", factory=SimpleDB)
     dbAe.meta['tip'] = "Spatial database for power-law coefficient, Ae."
 
+    dbStateVars = pythia.pyre.inventory.facility("db_state_variables", family="spatial_database", factory=SimpleDB)
+    dbStateVars.meta["tip"] = "Spatial database for power-law state variables, viscous_strain, reference_stress, reference_strain"
+
     from spatialdata.spatialdb.generator.Geometry import Geometry
     geometry = pythia.pyre.inventory.facility("geometry", family="geometry", factory=Geometry)
     geometry.meta['tip'] = "Geometry for output database."
@@ -81,71 +85,197 @@ class PowerLawApp(Application):
         coordsys = self.geometry.coordsys
 
         (npoints, spaceDim) = points.shape
-        refStrainRate = numpy.zeros((npoints,), dtype=numpy.float64)
-        refStress = numpy.zeros((npoints,), dtype=numpy.float64)
-
+        pwerlawRefStrainRate = numpy.zeros((npoints,), dtype=numpy.float64)
+        powerlawRefStress = numpy.zeros((npoints,), dtype=numpy.float64)
+        
         # Query databases to get inputs at output points
         self._info.log("Querying for parameters at output points.")
-        n = self._queryDB(self.dbExponent, "power-law-exponent", points, coordsys)
-        Q = self._queryDB(self.dbActivationE, "activation-energy", points, coordsys)
-        logAe = self._queryDB(self.dbAe, "log-flow-constant", points, coordsys)
-        scaleAe = self._queryDB(self.dbAe, "flow-constant-scale", points, coordsys)
-        T = self._queryDB(self.dbTemperature, "temperature", points, coordsys)
+        n = self._queryDB(self.dbExponent, ["power_law_exponent"], points, coordsys)
+        Q = self._queryDB(self.dbActivationE, ["activation_energy"], points, coordsys)
+        logAe = self._queryDB(self.dbAe, ["log_flow_constant"], points, coordsys)
+        scaleAe = self._queryDB(self.dbAe, ["flow_constant_scale"], points, coordsys)
+        T = self._queryDB(self.dbTemperature, ["temperature"], points, coordsys)
 
         # Compute power-law parameters
-        self._info.log("Computing parameters at output points.")
+        self._info.log("Computing parameters at output points...")
         from pythia.pyre.handbook.constants.fundamental import R
         Ae = 10**(logAe - scaleAe * n)
         At = 3**(0.5 * (n + 1)) / 2.0 * Ae * numpy.exp(-Q / (R.value * T))
 
         if self.refSelection == "stress":
-            refStress[:] = self.refStress.value
-            refStrainRate = self.refStress.value**n * At
+            powerlawRefStress[:] = self.powerlawRefStress.value
+            pwerlawRefStrainRate = self.powerlawRefStress.value**n * At
         elif self.refSelection == "strain_rate":
-            refStrainRate[:] = self.refStrainRate.value
-            refStress = (self.refStrainRate.value / At)**(1.0 / n)
+            pwerlawRefStrainRate[:] = self.pwerlawRefStrainRate.value
+            powerlawRefStress = (self.pwerlawRefStrainRate.value / At)**(1.0 / n)
         else:
             raise ValueError(f"Invalid value (self.refSelection) for reference value.")
 
-        refStressInfo = {
-            'name': "reference-stress",
+        powerlawRefStressInfo = {
+            'name': "power_law_reference_stress",
             'units': "Pa",
-            'data': refStress.flatten()
+            'data': powerlawRefStress.flatten()
         }
-        refStrainRateInfo = {
-            'name': "reference-strain-rate",
+        pwerlawRefStrainRateInfo = {
+            'name': "power_law_reference_strain_rate",
             'units': "1/s",
-            'data': refStrainRate.flatten()
+            'data': pwerlawRefStrainRate.flatten()
         }
         exponentInfo = {
-            'name': "power-law-exponent",
+            'name': "power_law_exponent",
             'units': "none",
             'data': n.flatten()
         }
 
+        tensorComponents = self.TENSOR_COMPONETS if spaceDim == 3 else self.TENSOR_COMPONENTS[:4]
+
+        valueNames = [f"viscous_strain_{component}" for component in tensorComponents]
+        viscousStrain = self._queryDB(self.dbStateVars, valueNames, points, coordsys)
+        viscousStrainInfo = [{
+            "name": "viscous_strain_xx",
+            "units": "none",
+            "data": viscousStrain[:,0].flatten(),
+            },{
+            "name": "viscous_strain_yy",
+            "units": "none",
+            "data": viscousStrain[:,1].flatten(),
+            },{
+            "name": "viscous_strain_zz",
+            "units": "none",
+            "data": viscousStrain[:,2].flatten(),
+            },{
+            "name": "viscous_strain_xy",
+            "units": "none",
+            "data": viscousStrain[:,3].flatten(),
+            }]
+        if spaceDim == 3:
+            viscousStrainInfo += [{
+                "name": "viscous_strain_yz",
+                "units": "none",
+                "data": viscousStrain[:,4].flatten(),
+                },{
+                "name": "viscous_strain_xz",
+                "units": "none",
+                "data": viscousStrain[:,5].flatten(),
+                }]
+
+        valueNames = [f"deviatoric_stress_{component}" for component in tensorComponents]
+        devStress = self._queryDB(self.dbStateVars, valueNames, points, coordsys)
+        devStressInfo = [{
+            "name": "deviatoric_stress_xx",
+            "units": "none",
+            "data": viscousStrain[:,0].flatten(),
+            },{
+            "name": "deviatoric_stress_yy",
+            "units": "none",
+            "data": viscousStrain[:,1].flatten(),
+            },{
+            "name": "deviatoric_stress_zz",
+            "units": "none",
+            "data": viscousStrain[:,2].flatten(),
+            },{
+            "name": "deviatoric_stress_xy",
+            "units": "none",
+            "data": viscousStrain[:,3].flatten(),
+            }]
+        if spaceDim == 3:
+            viscousStrainInfo += [{
+                "name": "deviatoric_stress_yz",
+                "units": "none",
+                "data": viscousStrain[:,4].flatten(),
+                },{
+                "name": "deviatoric_stress_xz",
+                "units": "none",
+                "data": viscousStrain[:,5].flatten(),
+                }]
+
+        valueNames = [f"reference_stress_{component}" for component in tensorComponents]
+        refStress = self._queryDB(self.dbStateVars, valueNames, points, coordsys)
+        refStressInfo = [{
+            "name": "reference_stress_xx",
+            "units": "Pa",
+            "data": refStress[:,0].flatten(),
+            },{
+            "name": "reference_stress_yy",
+            "units": "Pa",
+            "data": refStress[:,1].flatten(),
+            },{
+            "name": "reference_stress_zz",
+            "units": "Pa",
+            "data": refStress[:,2].flatten(),
+            },{
+            "name": "reference_stress_xy",
+            "units": "Pa",
+            "data": refStress[:,3].flatten(),
+            }]
+        if spaceDim == 3:
+            refStressInfo += [{
+                "name": "reference_stress_yz",
+                "units": "Pa",
+                "data": refStress[:,4].flatten(),
+                },{
+                "name": "reference_stress_xz",
+                "units": "Pa",
+                "data": refStress[:,5].flatten(),
+                }]
+
+        valueNames = [f"reference_strain_{component}" for component in tensorComponents]
+        refStrain = self._queryDB(self.dbStateVars, valueNames, points, coordsys)
+        refStrainInfo = [{
+            "name": "reference_strain_xx",
+            "units": "none",
+            "data": refStrain[:,0].flatten(),
+            },{
+            "name": "reference_strain_yy",
+            "units": "none",
+            "data": refStrain[:,1].flatten(),
+            },{
+            "name": "reference_strain_zz",
+            "units": "none",
+            "data": refStrain[:,2].flatten(),
+            },{
+            "name": "reference_strain_xy",
+            "units": "none",
+            "data": refStrain[:,3].flatten(),
+            }]
+        if spaceDim == 3:
+            refStrainInfo += [{
+                "name": "reference_strain_yz",
+                "units": "none",
+                "data": refStrain[:,4].flatten(),
+                },{
+                "name": "reference_strain_xz",
+                "units": "none",
+                "data": refStrain[:,5].flatten(),
+                }]
+
         # Write database
         self._info.log("Writing database.")
         data = {
             'points': points,
             'coordsys': coordsys,
             'data_dim': self.geometry.dataDim,
-            'values': [refStressInfo, refStrainRateInfo, exponentInfo]
+            'values': [powerlawRefStressInfo, pwerlawRefStrainRateInfo, exponentInfo]
         }
+        data["values"] += viscousStrainInfo
+        data["values"] += refStressInfo
+        data["values"] += refStrainInfo
+        data["values"] += devStressInfo
         from spatialdata.spatialdb.SimpleIOAscii import createWriter
         writer = createWriter(self.dbFilename)
         writer.write(data)
 
-    def _queryDB(self, db, valueName, points, cs):
+    def _queryDB(self, db, valueNames, points, cs):
         """
         Query spatial database
         """
 
         (npoints, spaceDim) = points.shape
-        data = numpy.zeros((npoints, 1), dtype=numpy.float64)
-        err = numpy.zeros((npoints,), dtype=numpy.int64)
+        data = numpy.zeros((npoints, len(valueNames)), dtype=numpy.float64)
+        err = numpy.zeros((npoints,), dtype=numpy.int32)
 
         db.open()
-        db.setQueryValues([valueName])
+        db.setQueryValues(valueNames)
         db.multiquery(data, err, points, cs)
         db.close()
         errSum = numpy.sum(err)

docs/user/examples/reverse-2d/step08-twofaults-powerlaw.md

@@ -1,3 +1,4 @@
+(sec-user-examples-reverse-2d-step08)=
 # Step 8: Slip on Two Faults and Power-law Viscoelastic Materials
 
 % Metadata extracted from parameter files.
@@ -18,13 +19,44 @@ caption: Power-law viscoelastic bulk rheology parameters for Step 8.
 slab.bulk_rheology = pylith.materials.IsotropicPowerLaw
 
 [pylithapp.problem.materials.slab]
-db_auxiliary_field = spatialdata.spatialdb.SimpleDB
-db_auxiliary_field.description = Maxwell viscoelastic properties
-db_auxiliary_field.iohandler.filename = mat_powerlaw.spatialdb
+db_auxiliary_field = spatialdata.spatialdb.CompositeDB
+db_auxiliary_field.description = Power law material properties
 
 bulk_rheology.auxiliary_subfields.power_law_reference_strain_rate.basis_order = 0
 bulk_rheology.auxiliary_subfields.power_law_reference_stress.basis_order = 0
 bulk_rheology.auxiliary_subfields.power_law_exponent.basis_order = 0
+
+[pylithapp.problem.materials.slab.db_auxiliary_field]
+# Elastic properties
+values_A = [density, vs, vp]
+db_A = spatialdata.spatialdb.SimpleDB
+db_A.description = Elastic properties for slab
+db_A.iohandler.filename = mat_elastic.spatialdb
+
+# Power law properties
+values_B = [
+	 power_law_reference_stress,  power_law_reference_strain_rate,  power_law_exponent,
+	 viscous_strain_xx, viscous_strain_yy, viscous_strain_zz, viscous_strain_xy,
+	 reference_stress_xx, reference_stress_yy, reference_stress_zz, reference_stress_xy,
+	 reference_strain_xx, reference_strain_yy, reference_strain_zz, reference_strain_xy,
+	 deviatoric_stress_xx,  deviatoric_stress_yy,  deviatoric_stress_zz,  deviatoric_stress_xy
+	 ]
+db_B = spatialdata.spatialdb.SimpleDB
+db_B.description = Material properties specific to power law bulk rheology for the slab
+db_B.iohandler.filename = mat_powerlaw.spatialdb
+db_B.query_type = linear
+```
+
+## Power-law spatial database
+
+We use the utility script `pylith_powerlaw_gendb` (see {ref}`sec-user-run-pylith-pylith-powerlaw-gendb`) to generate the spatial database `mat_powerlaw.spatialdb` with the power-law bulk rheology parameters.
+We provide the parameters for `pylith_powerlaw_gendb` in `powerlaw_gendb.cfg`, which follows the same formatting conventions as the PyLith parameter files. 
+
+```{code-block} console
+---
+caption: Generate spatial database with power-law viscoelastic material properties.
+---
+$ pylith_powerlaw_gendb powerlaw_gendb.cfg
 ```
 
 ## Running the simulation
@@ -48,11 +80,11 @@ $ pylith step08_twofaults_powerlaw.cfg
 # -- many lines omitted --
 
 25 TS dt 0.2 time 4.8
-    0 SNES Function norm 5.602869611772e-06 
-    Linear solve converged due to CONVERGED_ATOL iterations 262
-    1 SNES Function norm 6.656502817834e-08 
-    Linear solve converged due to CONVERGED_ATOL iterations 117
-    2 SNES Function norm 9.645013365143e-10 
+    0 SNES Function norm 2.142894498538e-06
+    Linear solve converged due to CONVERGED_ATOL iterations 200
+    1 SNES Function norm 6.320401896875e-09
+    Linear solve converged due to CONVERGED_ATOL iterations 67
+    2 SNES Function norm 1.048498421861e-10
   Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 2
 26 TS dt 0.2 time 5.
  >> /software/baagaard/py38-venv/pylith-debug/lib/python3.8/site-packages/pylith/problems/Problem.py:201:finalize

docs/user/examples/strikeslip-2d/step04-varslip.md

@@ -43,7 +43,7 @@ defaults.quadrature_order = 2
 
 [pylithapp.problem.solution.subfields]
 displacement.basis_order = 2
-lagrange_fault.basis_order = 2
+lagrange_multiplier_fault.basis_order = 2
 
 [pylithapp.problem]
 solution_observers = [domain, top_boundary, bot_boundary, gps_stations]

docs/user/run-pylith/utilities.md

@@ -299,9 +299,6 @@ An additional parameter defines the units of the activation energy.
 You must also specify either a reference stress or a reference strain rate.
 
 You place all of the application parameters in `powerlaw_gendb.cfg`, which the application will read by default.
-
-:::{admonition} TODO
-Add cross reference to example using `pylith_powerlaw_gendb`.
-:::
+See {ref}`sec-user-examples-reverse-2d-step08` for an example of how to use `pylith_powerlaw_gendb`.
 
 % End of file

examples/barwaves-2d/step04_swave_prescribedslip.cfg

@@ -62,7 +62,7 @@ defaults.quadrature_order = 1
 [pylithapp.problem.solution.subfields]
 displacement.basis_order = 1
 velocity.basis_order = 1
-lagrange_fault.basis_order = 1
+lagrange_multiplier_fault.basis_order = 1
 
 # ----------------------------------------------------------------------
 # fault

examples/meshing-cubit/surface-nurbs/merge-surfs/pylithapp.cfg

@@ -42,8 +42,8 @@ solution = pylith.problems.SolnDispLagrange
 displacement.basis_order = 1
 displacement.quadrature_order = 1
 
-lagrange_fault.basis_order = 1
-lagrange_fault.quadrature_order = 1
+lagrange_multiplier_fault.basis_order = 1
+lagrange_multiplier_fault.quadrature_order = 1
 
 [pylithapp.problem]
 solution_observers = [domain]

examples/reverse-2d/Makefile.am

@@ -43,6 +43,10 @@ dist_noinst_DATA = \
 	mat_gravity_refstate.spatialdb \
 	mat_maxwell.spatialdb \
 	mat_powerlaw.spatialdb \
+	powerlaw_params.spatialdb \
+	powerlaw_temperature.spatialdb \
+	powerlaw_gendb.cfg \
+	powerlaw_points.txt \
 	traction_surfload.spatialdb \
 	viz/plot_dispwarp.py \
 	viz/plot_mesh.py