Merge pull request #741 from baagaard-usgs/fix-examples

Small fixes to strikeslip-2d and reverse-2d examples

docs/user/examples/reverse-2d/step01-gravity.md

@@ -312,7 +312,7 @@ Increasing the basis order reduces the maximum shear stress by about 10 orders o
 
 1. The displacement field is very similar for the three different discretizations.
 2. In this case, we know the vertical displacement in the exact solution depends on the square of the depth, so we expect the numerical accuracy to be limited when we use a basis order of 1 for the displacement solution subfield.
-3. In the exact solution, the axial components of the Cauchy stress tensor vary linear with depth, and the xy component of the Cauchy stress is zero.
+3. In the exact solution, the axial components of the Cauchy stress tensor increase linearly with depth, and the xy component of the Cauchy stress is zero.
 4. With the coarse mesh, the shear stress is close to 10 MPa throughout the mesh.
 5. Refining the mesh decreases the magnitude shear stress by about a factor of 2.
 6. Using a basis order of 2 for the solution subfields reduces the magnitude of the shear stress by nearly 10 orders of magnitude.

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

@@ -24,6 +24,29 @@ Boundary conditions for static coseismic slip.
 We set the x and y displacement to zero on the +x and -x boundaries and prescribe left-lateral slip that varies along strike.
 :::
 
+For greater accuracy in modeling the spatial variation in slip, we use a basis order of 2 for the solution subfields.
+
+```{code-block} cfg
+---
+caption: Parameters related to increasing the basis order of the solution subfields to 2.
+---
+[pylithapp.problem]
+defaults.quadrature_order = 2
+
+[pylithapp.problem.solution.subfields]
+displacement.basis_order = 2
+lagrange_multiplier_fault.basis_order = 2
+
+
+[pylithapp.problem.materials.elastic_xneg]
+derived_subfields.cauchy_strain.basis_order = 1
+derived_subfields.cauchy_stress.basis_order = 1
+
+[pylithapp.problem.materials.elastic_xpos]
+derived_subfields.cauchy_strain.basis_order = 1
+derived_subfields.cauchy_stress.basis_order = 1
+```
+
 We also add output of the solution at fake GNSS stations given in the file `gnss_stations.txt`.
 You can use the Python script `generate_gnssstations.py` to generate a different random set of stations; the default parameters will generate the provided `gnss_stations.txt` file.
 
@@ -47,6 +70,20 @@ reader.filename = gnss_stations.txt
 reader.coordsys.space_dim = 2
 ```
 
+```{code-block} cfg
+---
+caption: Solution and output parameters for Step 4. We add output of the solution at fake GNSS stations.
+---
+[pylithapp.problem]
+solution_observers = [domain, top_boundary, bot_boundary, gnss_stations]
+solution_observers.gnss_stations = pylith.meshio.OutputSolnPoints
+
+[pylithapp.problem.solution_observers.gnss_stations]
+label = gnss_stations
+reader.filename = gnss_stations.txt
+reader.coordsys.space_dim = 2
+```
+
 The earthquake rupture occurs along the central portion of the fault with spatially variable slip.
 
 :::{figure-md} fig:example:strikeslip:2d:step04:slip
@@ -96,10 +133,10 @@ $ pylith step04_varslip.cfg
  -- timedependent(info)
  -- Solving problem.
 0 TS dt 0.01 time 0.
-    0 SNES Function norm 3.786132692381e-03
-    Linear solve converged due to CONVERGED_ATOL iterations 61
-    1 SNES Function norm 2.070681594908e-12
-  Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
+    0 SNES Function norm 5.220560316093e-03
+      Linear solve converged due to CONVERGED_ATOL iterations 27
+    1 SNES Function norm 1.523809186100e-12
+    Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
 1 TS dt 0.01 time 0.01
  >> /software/unix/py3.12-venv/pylith-debug/lib/python3.12/site-packages/pylith/problems/Problem.py:199:finalize
  -- timedependent(info)
@@ -108,7 +145,7 @@ $ pylith step04_varslip.cfg
 
 The beginning of the output written to the terminal matches that in our previous simulations.
 At the end of the output written to the terminal, we see that the solver advanced the solution one time step (static simulation).
-The linear solve converged after 61 iterations and the norm of the residual met the absolute convergence tolerance (`ksp_atol`).
+The linear solve converged after 27 iterations and the norm of the residual met the absolute convergence tolerance (`ksp_atol`).
 The nonlinear solve converged in 1 iteration, which we expect because this is a linear problem, and the residual met the absolute convergence tolerance (`snes_atol`).
 
 ## Visualizing the results

examples/strikeslip-2d/pylithapp.cfg

@@ -164,14 +164,14 @@ label = boundary_xpos
 label_value = 11
 constrained_dof = [0, 1]
 db_auxiliary_field = pylith.bc.ZeroDB
-db_auxiliary_field.description = Dirichlet BC +x boundary
+db_auxiliary_field.description = Dirichlet BC on +x boundary
 
 [pylithapp.problem.bc.bc_xneg]
 label = boundary_xneg
 label_value = 10
 constrained_dof = [0, 1]
 db_auxiliary_field = pylith.bc.ZeroDB
-db_auxiliary_field.description = Dirichlet BC -x boundary
+db_auxiliary_field.description = Dirichlet BC on -x boundary
 
 
 # End of file

examples/strikeslip-2d/step04_varslip.cfg

@@ -22,7 +22,7 @@
 base = [pylithapp.cfg]
 description = Coseismic prescribed slip with zero displacement Dirichlet boundary conditions.
 authors = [Brad Aagaard]
-keywords = [prescribed slip]
+keywords = [prescribed slip, basis order 2]
 arguments = [step04_varslip.cfg]
 version = 1.0.0
 pylith_version = [>=4.0, <5.0]
@@ -44,15 +44,33 @@ problem.progress_monitor.filename = output/step04_varslip-progress.txt
 # output filenames. The default directory for output is 'output'.
 problem.defaults.name = step04_varslip
 
-# ----------------------------------------------------------------------
-# mesh_generator
-# ----------------------------------------------------------------------
-[pylithapp.mesh_generator]
-refiner = pylith.topology.RefineUniform
-
 # ----------------------------------------------------------------------
 # problem
 # ----------------------------------------------------------------------
+# Update the parameters to use a basis order of 2 for the solution fields.
+#
+# The accuracy of the stress and strain will be 1 order lower than the basis
+# order of the displacement field. Consequently, we use a basis order of 1
+# (rather than 0) for the output of the Cauchy stress and strain.
+#
+# We do not change the basis order of the output of the displacement field,
+# because many visualization tools do not know how to display fields with
+# a basis order of 2.
+
+[pylithapp.problem]
+# Set the default quadrature order for all subfields. We want to set it
+# to the maximum of the basis order of the solution subfields.
+#
+# IMPORTANT: The quadrature order *must* the same for all solution and
+# auxiliary subfields. PyLith will verify that this requirement is met.
+# This requirement may be relaxed in the future.
+defaults.quadrature_order = 2
+
+[pylithapp.problem.solution.subfields]
+displacement.basis_order = 2
+lagrange_multiplier_fault.basis_order = 2
+
+
 [pylithapp.problem]
 # We add output at our fake GNSS stations that we will use a fake observations.
 solution_observers = [domain, top_boundary, bot_boundary, gnss_stations]
@@ -65,6 +83,18 @@ reader.filename = gnss_stations.txt
 reader.coordsys.space_dim = 2
 
 
+# ----------------------------------------------------------------------
+# materials
+# ----------------------------------------------------------------------
+[pylithapp.problem.materials.elastic_xneg]
+derived_subfields.cauchy_strain.basis_order = 1
+derived_subfields.cauchy_stress.basis_order = 1
+
+[pylithapp.problem.materials.elastic_xpos]
+derived_subfields.cauchy_strain.basis_order = 1
+derived_subfields.cauchy_stress.basis_order = 1
+
+
 # ----------------------------------------------------------------------
 # fault
 # ----------------------------------------------------------------------

examples/strikeslip-2d/viz/plot_inversion_results.py

@@ -20,7 +20,7 @@ class PlotApp:
 
 
     COLOR_OBSERVED = "black"
-    COLORS_MODEL = ["blue", "orange", "purple", "red", "cyan"]
+    COLORS_MODEL = ["orange", "purple", "red", "cyan"]
     XLIM = (-25.0, 25.0)
 
     def main(self, filename_models: str=None, filename_observed: str=None, observed_only: bool=False, show_plot: bool=False):
@@ -78,7 +78,7 @@ def _plot_models(self, axes):
         nmodels = self.models_slip.shape[1]
         coords = self.models_coords / 1.0e+3
         for imodel in range(nmodels):
-            axes.plot(coords, self.models_slip[:,imodel], color=self.COLORS_MODEL[imodel], label=f"Model penalty={self.models_penalty[imodel]}")
+            axes.plot(coords, self.models_slip[:,imodel], "--", color=self.COLORS_MODEL[imodel], label=f"Model penalty={self.models_penalty[imodel]}")
 
 
 def cli():

release-notes/checklist.txt

@@ -37,6 +37,7 @@ DISTRIBUTIONS
   * Generate binaries using pylith_installer packager/build.py.
     Check on various platforms.
     Check trapping of errors.
+    Check macOSX minos `for i in lib*.dylib; do otool -l $i | grep -E -A4 '(LC_VERSION_MIN_MACOSX|LC_BUILD_VERSION)' | grep -B1 sdk; done`
 
 TAG