diff --git a/CHANGES.md b/CHANGES.md index 695a76079b..c5b292dbfa 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -3,8 +3,15 @@ See 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 diff --git a/applications/pylith_powerlaw_gendb b/applications/pylith_powerlaw_gendb index a3633d1286..77f9b351d4 100755 --- a/applications/pylith_powerlaw_gendb +++ b/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) diff --git a/docs/user/examples/reverse-2d/step08-twofaults-powerlaw.md b/docs/user/examples/reverse-2d/step08-twofaults-powerlaw.md index 942ade5866..1ca0aef807 100644 --- a/docs/user/examples/reverse-2d/step08-twofaults-powerlaw.md +++ b/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 diff --git a/docs/user/examples/strikeslip-2d/step04-varslip.md b/docs/user/examples/strikeslip-2d/step04-varslip.md index 605373af3f..268633a687 100644 --- a/docs/user/examples/strikeslip-2d/step04-varslip.md +++ b/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] diff --git a/docs/user/run-pylith/utilities.md b/docs/user/run-pylith/utilities.md index 7cb0d0fbb6..3e0b799deb 100644 --- a/docs/user/run-pylith/utilities.md +++ b/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 diff --git a/examples/barwaves-2d/step04_swave_prescribedslip.cfg b/examples/barwaves-2d/step04_swave_prescribedslip.cfg index 153a806eff..6f618cf6ce 100644 --- a/examples/barwaves-2d/step04_swave_prescribedslip.cfg +++ b/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 diff --git a/examples/meshing-cubit/surface-nurbs/merge-surfs/pylithapp.cfg b/examples/meshing-cubit/surface-nurbs/merge-surfs/pylithapp.cfg index 328de93150..766382c9bf 100644 --- a/examples/meshing-cubit/surface-nurbs/merge-surfs/pylithapp.cfg +++ b/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] diff --git a/examples/reverse-2d/Makefile.am b/examples/reverse-2d/Makefile.am index f6f22bc12d..90de5f6b3a 100644 --- a/examples/reverse-2d/Makefile.am +++ b/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 diff --git a/examples/reverse-2d/mat_powerlaw.spatialdb b/examples/reverse-2d/mat_powerlaw.spatialdb index b5687d9459..c1e8524683 100644 --- a/examples/reverse-2d/mat_powerlaw.spatialdb +++ b/examples/reverse-2d/mat_powerlaw.spatialdb @@ -1,39 +1,44 @@ #SPATIAL.ascii 1 SimpleDB { - num-values = 22 - value-names = density vs vp power_law_reference_strain_rate power_law_reference_stress power_law_exponent viscous_strain_xx viscous_strain_yy viscous_strain_zz viscous_strain_xy deviatoric_stress_xx deviatoric_stress_yy deviatoric_stress_zz deviatoric_stress_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 - value-units = kg/m**3 m/s m/s 1/s Pa None None None None None Pa Pa Pa Pa Pa Pa Pa Pa None None None None - num-locs = 1 - data-dim = 0 + num-values = 19 + value-names = 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 + value-units = Pa 1/s none none none none none Pa Pa Pa Pa none none none none none none none none + num-locs = 31 + data-dim = 1 space-dim = 2 cs-data = cartesian { - to-meters = 1 - space-dim = 2 - } + to-meters = 1 + space-dim = 2 } -// Columns are -// (1) x coordinate (m) -// (2) y coordinate (m) -// (3) density (kg/m**3) -// (4) vs (m/s) -// (5) vp (m/s) -// (6) Power-law reference stress (Pa) -// (7) Power-law reference strain rate (1/s) -// (8) Power-law exponent (None) -// (9) viscous_strain_xx (None) -// (10) viscous_strain_yy (None) -// (11) viscous_strain_zz (None) -// (12) viscous_strain_xy (None) -// (13) deviatoric_stress_xx (Pa) -// (14) deviatoric_stress_yy (Pa) -// (15) deviatoric_stress_zz (Pa) -// (16) deviatoric_stress_xy (Pa) -// (17) reference_stress_xx (Pa) -// (18) reference_stress_yy (Pa) -// (19) reference_stress_zz (Pa) -// (29) reference_stress_xy (Pa) -// (21) reference_strain_xx (None) -// (22) reference_strain_yy (None) -// (23) reference_strain_zz (None) -// (24) reference_strain_xy (None) - 0.0 0.0 2500.0 3464.1016 6000.0 1.000000e-06 7.0e+8 3.500 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 +} + 0.000000e+00 0.000000e+00 1.818610e+16 1.000000e-06 1.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.000000e+02 7.635317e+14 1.000000e-06 1.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.000000e+03 5.936971e+13 1.000000e-06 1.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.500000e+03 7.258145e+12 1.000000e-06 1.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.900000e+03 1.735232e+12 1.000000e-06 1.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.100000e+03 6.971431e+16 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.500000e+03 9.728762e+15 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.000000e+03 1.131492e+15 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.500000e+03 1.743727e+14 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.000000e+03 3.380856e+13 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.500000e+03 7.924871e+12 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.000000e+03 2.176987e+12 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.500000e+03 6.837569e+11 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.000000e+03 2.407282e+11 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.500000e+03 9.348560e+10 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.000000e+03 3.952049e+10 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.500000e+03 1.798919e+10 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.000000e+03 8.736679e+09 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.500000e+03 4.492496e+09 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.000000e+03 2.429991e+09 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.500000e+03 1.374933e+09 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -2.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -3.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -4.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -5.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -6.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -7.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -8.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -9.000000e+04 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 + 0.000000e+00 -1.000000e+05 8.099233e+08 1.000000e-06 3.500000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 diff --git a/examples/reverse-2d/powerlaw_gendb.cfg b/examples/reverse-2d/powerlaw_gendb.cfg new file mode 100644 index 0000000000..89877a5979 --- /dev/null +++ b/examples/reverse-2d/powerlaw_gendb.cfg @@ -0,0 +1,70 @@ +# Parameter file for generating spatial database for power law viscoelastic +# bulk rheology using `pylith_powerlaw_gendb`. +[powerlaw_gendb] +# Output filename +database_filename = mat_powerlaw.spatialdb + +# Use strain rate as the reference for computing power law parameters. +reference_value = strain_rate +reference_strain_rate = 1.0e-6/s + +db_state_variables = spatialdata.spatialdb.UniformDB + + +# Spatial databaseses for power law parameters. +[powerlaw_gendb.db_exponent] +description = Power-law exponent +iohandler.filename = powerlaw_params.spatialdb + +[powerlaw_gendb.db_activation_energy] +description = Activation energy +iohandler.filename = powerlaw_params.spatialdb + +[powerlaw_gendb.db_powerlaw_coefficient] +description = Experimentally derived power-law coefficient +iohandler.filename = powerlaw_params.spatialdb + +[powerlaw_gendb.db_temperature] +# Spatial database for temperature profile. +description = Temperature +query_type = linear +iohandler.filename = powerlaw_temperature.spatialdb + +[powerlaw_gendb.db_state_variables] +# Spatial database for reference stress and strain state and state variables. +description = Power-law state variables +values = [ + viscous_strain_xx, + viscous_strain_yy, + viscous_strain_zz, + viscous_strain_xy, + deviatoric_stress_xx, + deviatoric_stress_yy, + deviatoric_stress_zz, + deviatoric_stress_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, + ] +data = [ + 0.0, 0.0, 0.0, 0.0, + 0.0*Pa, 0.0*Pa, 0.0*Pa, 0.0*Pa, + 0.0*Pa, 0.0*Pa, 0.0*Pa, 0.0*Pa, + 0.0, 0.0, 0.0, 0.0, + ] + +[powerlaw_gendb.geometry] +# 1D vertical profile of points +coordsys.space_dim = 2 +data_dim = 1 + +reader = spatialdata.utils.PointsStream +reader.filename = powerlaw_points.txt + + +# End of file \ No newline at end of file diff --git a/examples/reverse-2d/powerlaw_params.spatialdb b/examples/reverse-2d/powerlaw_params.spatialdb new file mode 100644 index 0000000000..9dbd57b15c --- /dev/null +++ b/examples/reverse-2d/powerlaw_params.spatialdb @@ -0,0 +1,37 @@ +// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting) +// +// This spatial database defines the power-law parameters for two materials: +// Wet granite and dry olivine. The properties are from Strehlau and +// Meissner (1987). Wet granite is specified from zero to 2 km depth, +// while dry olivine is specified from 2 to 100 km depth. +// Note that the flow constant is expressed as log10(flow-constant), where +// the units are GPa^(-n)/s. The activation energy (Q) is expressed in +// kJ/mol. To accommodate the strange units of the flow constant, the +// given units of are 'None', and we specify a new parameter +// (flow-constant-multiplier), with a value of 9, to indicate that the +// underlying units are GPa. A value of 6 would indicate MPa, etc. +// +#SPATIAL.ascii 1 +SimpleDB { + num-values = 4 // number of material property values + value-names = log_flow_constant activation_energy power_law_exponent flow_constant_scale // names of the material property values + value-units = None kJ/mol None None // units + num-locs = 4 // number of locations + data-dim = 1 + space-dim = 2 + cs-data = cartesian { + to-meters = 1000.0 + space-dim = 2 + } +} +// Columns are +// (1) x coordinate (m) +// (2) y coordinate (m) +// (3) flow_constant (None) +// (4) activation_energy (kJ/mol) +// (5) power_law_exponent (None) +// (6) activation_energy_multiplier (None) +0.0 0.0 2.0 137.0 1.5 9.0 +0.0 -1.99 2.0 137.0 1.5 9.0 +0.0 -2.01 15.5 535.0 3.5 9.0 +0.0 -100.0 15.5 535.0 3.5 9.0 diff --git a/examples/reverse-2d/powerlaw_points.txt b/examples/reverse-2d/powerlaw_points.txt new file mode 100644 index 0000000000..fc43e42a26 --- /dev/null +++ b/examples/reverse-2d/powerlaw_points.txt @@ -0,0 +1,32 @@ +# Locations at which we want to output power-law properties. +0.0 0.0 +0.0 -500.0 +0.0 -1000.0 +0.0 -1500.0 +0.0 -1900.0 +0.0 -2100.0 +0.0 -2500.0 +0.0 -3000.0 +0.0 -3500.0 +0.0 -4000.0 +0.0 -4500.0 +0.0 -5000.0 +0.0 -5500.0 +0.0 -6000.0 +0.0 -6500.0 +0.0 -7000.0 +0.0 -7500.0 +0.0 -8000.0 +0.0 -8500.0 +0.0 -9000.0 +0.0 -9500.0 +0.0 -10000.0 +0.0 -20000.0 +0.0 -30000.0 +0.0 -40000.0 +0.0 -50000.0 +0.0 -60000.0 +0.0 -70000.0 +0.0 -80000.0 +0.0 -90000.0 +0.0 -100000.0 diff --git a/examples/reverse-2d/powerlaw_temperature.spatialdb b/examples/reverse-2d/powerlaw_temperature.spatialdb new file mode 100644 index 0000000000..492f3c8e07 --- /dev/null +++ b/examples/reverse-2d/powerlaw_temperature.spatialdb @@ -0,0 +1,26 @@ +// -*- C++ -*- (tell Emacs to use C++ mode for syntax highlighting) +// +// This spatial database defines a simple picewise linear 1D temperature profile, +// assuming a temperature of 100 degrees C (373 degrees K) at the surface, +// increasing to 1,000 degrees C (1273 degrees K) at 10 km depth. +// +#SPATIAL.ascii 1 +SimpleDB { + num-values = 1 // number of material property values + value-names = temperature // names of the material property values + value-units = K // units + num-locs = 3 // number of locations + data-dim = 1 + space-dim = 2 + cs-data = cartesian { + to-meters = 1000.0 + space-dim = 2 + } +} +// Columns are +// (1) x coordinate (m) +// (2) y coordinate (m) +// (3) temperature (K) +0.0 0.0 373.0 +0.0 -10.0 1273.0 +0.0 -100.0 1273.0 diff --git a/examples/reverse-2d/pylithapp.cfg b/examples/reverse-2d/pylithapp.cfg index 416ab1b05a..87dbede439 100644 --- a/examples/reverse-2d/pylithapp.cfg +++ b/examples/reverse-2d/pylithapp.cfg @@ -72,7 +72,6 @@ label_value = 1 db_auxiliary_field = spatialdata.spatialdb.SimpleDB db_auxiliary_field.description = Elastic properties for slab -db_auxiliary_field.iohandler.filename = mat_elastic.spatialdb auxiliary_subfields.density.basis_order = 0 auxiliary_subfields.gravitational_acceleration.basis_order = 0 diff --git a/examples/reverse-2d/step01_gravity.cfg b/examples/reverse-2d/step01_gravity.cfg index 3fa3c3de85..145138446b 100644 --- a/examples/reverse-2d/step01_gravity.cfg +++ b/examples/reverse-2d/step01_gravity.cfg @@ -67,6 +67,8 @@ defaults.quadrature_order = 2 basis_order = 2 [pylithapp.problem.materials.slab] +db_auxiliary_field.iohandler.filename = mat_elastic.spatialdb + # We discretize the displacement field with a basis order of 2 # so the stress and strain computed from the displacement field # will have an accuracy of one order lower. diff --git a/examples/reverse-2d/step04_surfload.cfg b/examples/reverse-2d/step04_surfload.cfg index 56c9c4ee1b..c9ec169ac8 100644 --- a/examples/reverse-2d/step04_surfload.cfg +++ b/examples/reverse-2d/step04_surfload.cfg @@ -56,6 +56,12 @@ problem.progress_monitor.filename = output/step04_surfload-progress.txt # output filenames. The default directory for output is 'output'. problem.defaults.name = step04_surfload +# ---------------------------------------------------------------------- +# materials +# ---------------------------------------------------------------------- +[pylithapp.problem.materials.slab] +db_auxiliary_field.iohandler.filename = mat_elastic.spatialdb + # ---------------------------------------------------------------------- # boundary conditions # ---------------------------------------------------------------------- diff --git a/examples/reverse-2d/step05_onefault.cfg b/examples/reverse-2d/step05_onefault.cfg index d8713654aa..1ce6b970a5 100644 --- a/examples/reverse-2d/step05_onefault.cfg +++ b/examples/reverse-2d/step05_onefault.cfg @@ -55,6 +55,12 @@ problem.defaults.name = step05_onefault [pylithapp.problem] solution = pylith.problems.SolnDispLagrange +# ---------------------------------------------------------------------- +# materials +# ---------------------------------------------------------------------- +[pylithapp.problem.materials.slab] +db_auxiliary_field.iohandler.filename = mat_elastic.spatialdb + # ---------------------------------------------------------------------- # fault # ---------------------------------------------------------------------- diff --git a/examples/reverse-2d/step06_twofaults_elastic.cfg b/examples/reverse-2d/step06_twofaults_elastic.cfg index 5de5fe0fbb..eae9ca0e12 100644 --- a/examples/reverse-2d/step06_twofaults_elastic.cfg +++ b/examples/reverse-2d/step06_twofaults_elastic.cfg @@ -54,6 +54,12 @@ problem.defaults.name = step06_twofaults_elastic [pylithapp.problem] solution = pylith.problems.SolnDispLagrange +# ---------------------------------------------------------------------- +# materials +# ---------------------------------------------------------------------- +[pylithapp.problem.materials.slab] +db_auxiliary_field.iohandler.filename = mat_elastic.spatialdb + # ---------------------------------------------------------------------- # fault # ---------------------------------------------------------------------- diff --git a/examples/reverse-2d/step08_twofaults_powerlaw.cfg b/examples/reverse-2d/step08_twofaults_powerlaw.cfg index d88fae81a6..b7f7619d6b 100644 --- a/examples/reverse-2d/step08_twofaults_powerlaw.cfg +++ b/examples/reverse-2d/step08_twofaults_powerlaw.cfg @@ -32,7 +32,8 @@ features = [ pylith.materials.IsotropicPowerLaw, pylith.faults.FaultCohesiveKin, pylith.faults.KinSrcStep, - spatialdata.spatialdb.UniformDB + spatialdata.spatialdb.UniformDB, + spatialdata.spatialdb.CompositeDB ] # ---------------------------------------------------------------------- @@ -76,14 +77,36 @@ solution = pylith.problems.SolnDispLagrange 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 +# We use a composite spatial database with elastic properties from `mat_elastic.spatialdb` +# and the power-law properties from `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 + + # ---------------------------------------------------------------------- # fault # ---------------------------------------------------------------------- diff --git a/examples/strikeslip-2d/step04_varslip.cfg b/examples/strikeslip-2d/step04_varslip.cfg index 0fe101f969..e3d88ab7aa 100644 --- a/examples/strikeslip-2d/step04_varslip.cfg +++ b/examples/strikeslip-2d/step04_varslip.cfg @@ -57,7 +57,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] # We add output at our fake GPS stations that we will use a fake observations. diff --git a/examples/subduction-3d/step02_coseismic.cfg b/examples/subduction-3d/step02_coseismic.cfg index 6d9ab4e953..e65d2633ad 100644 --- a/examples/subduction-3d/step02_coseismic.cfg +++ b/examples/subduction-3d/step02_coseismic.cfg @@ -55,7 +55,7 @@ solution = pylith.problems.SolnDispLagrange [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 # ---------------------------------------------------------------------- # faults diff --git a/examples/subduction-3d/step03_interseismic.cfg b/examples/subduction-3d/step03_interseismic.cfg index 2060bd822b..288d452285 100644 --- a/examples/subduction-3d/step03_interseismic.cfg +++ b/examples/subduction-3d/step03_interseismic.cfg @@ -59,7 +59,7 @@ solution = pylith.problems.SolnDispLagrange [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 # ---------------------------------------------------------------------- # faults diff --git a/examples/subduction-3d/step04_eqcycle.cfg b/examples/subduction-3d/step04_eqcycle.cfg index 9703725ac0..c324ed8aa3 100644 --- a/examples/subduction-3d/step04_eqcycle.cfg +++ b/examples/subduction-3d/step04_eqcycle.cfg @@ -57,7 +57,7 @@ solution = pylith.problems.SolnDispLagrange [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 # ---------------------------------------------------------------------- # faults diff --git a/examples/subduction-3d/step06_slowslip.cfg b/examples/subduction-3d/step06_slowslip.cfg index 85d9aca96a..21e743dbe1 100644 --- a/examples/subduction-3d/step06_slowslip.cfg +++ b/examples/subduction-3d/step06_slowslip.cfg @@ -71,7 +71,7 @@ solution = pylith.problems.SolnDispLagrange [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 [pylithapp.problem] solution_observers = [domain, groundsurf, gps_stations] diff --git a/examples/subduction-3d/step07a_leftlateral.cfg b/examples/subduction-3d/step07a_leftlateral.cfg index ec8adac0ad..ef102c08f9 100644 --- a/examples/subduction-3d/step07a_leftlateral.cfg +++ b/examples/subduction-3d/step07a_leftlateral.cfg @@ -67,7 +67,7 @@ defaults.quadrature_order = 1 # for the solution fields should be the same as the basis order for the slip impulses; # higher basis orders are also valid. displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 [pylithapp.problem] # We add output at our fake GPS stations that we will use to invert for the slip. diff --git a/libsrc/pylith/fekernels/Poroelasticity.hh b/libsrc/pylith/fekernels/Poroelasticity.hh index 24512690aa..a13819d43e 100644 --- a/libsrc/pylith/fekernels/Poroelasticity.hh +++ b/libsrc/pylith/fekernels/Poroelasticity.hh @@ -618,35 +618,34 @@ public: #endif - //Calculate bulk density + // Calculate bulk density static inline void bulkDensity_asScalar(const PylithInt dim, - const PylithInt numS, - const PylithInt numA, - const PylithInt sOff[], - const PylithInt sOff_x[], - const PylithScalar s[], - const PylithScalar s_t[], - const PylithScalar s_x[], - const PylithInt aOff[], - const PylithInt aOff_x[], - const PylithScalar a[], - const PylithScalar a_t[], - const PylithScalar a_x[], - const PylithReal t, - const PylithScalar x[], - const PylithInt numConstants, - const PylithScalar constants[], - PylithReal* bulkDensity) { - + const PylithInt numS, + const PylithInt numA, + const PylithInt sOff[], + const PylithInt sOff_x[], + const PylithScalar s[], + const PylithScalar s_t[], + const PylithScalar s_x[], + const PylithInt aOff[], + const PylithInt aOff_x[], + const PylithScalar a[], + const PylithScalar a_t[], + const PylithScalar a_x[], + const PylithReal t, + const PylithScalar x[], + const PylithInt numConstants, + const PylithScalar constants[], + PylithReal* bulkDensity) { // Poroelastic Context pylith::fekernels::Poroelasticity::PoroelasticContext poroelasticContext; pylith::fekernels::Poroelasticity::setPoroelasticContextQS( &poroelasticContext, dim, numS, sOff, sOff_x, s, s_t, s_x, aOff, aOff_x, a, a_t, a_x, t, x); - bulkDensity = &(poroelasticContext.bulkDensity); + *bulkDensity = poroelasticContext.bulkDensity; - } //bulkDensity_asScalar + } // bulkDensity_asScalar // ============================================================================= // Velocity diff --git a/pylith/Makefile.am b/pylith/Makefile.am index 3ab011ec7d..bfe060b670 100644 --- a/pylith/Makefile.am +++ b/pylith/Makefile.am @@ -119,7 +119,9 @@ EXTRA_DIST = \ problems/SolnDispPres.py \ problems/SolnDispPresLagrange.py \ problems/SolnDispPresTracStrain.py \ + problems/SolnDispPresTracStrainLagrange.py \ problems/SolnDispPresTracStrainVelPdotTdot.py \ + problems/SolnDispPresTracStrainVelPdotTdotLagrange.py \ problems/SolnDispPresVel.py \ problems/SolnDispVel.py \ problems/SolnDispVelLagrange.py \ diff --git a/pylith/materials/AuxSubfieldsIsotropicLinearPoroelasticity.py b/pylith/materials/AuxSubfieldsIsotropicLinearPoroelasticity.py index 17de4b9deb..f5f4084362 100644 --- a/pylith/materials/AuxSubfieldsIsotropicLinearPoroelasticity.py +++ b/pylith/materials/AuxSubfieldsIsotropicLinearPoroelasticity.py @@ -53,6 +53,12 @@ class AuxSubfieldsIsotropicLinearPoroelasticity(PetscComponent): biotModulus = pythia.pyre.inventory.facility("biot_modulus", family="auxiliary_subfield", factory=Subfield) biotModulus.meta['tip'] = "Biot modulus subfield." + referenceStress = pythia.pyre.inventory.facility("reference_stress", family="auxiliary_subfield", factory=Subfield) + referenceStress.meta['tip'] = "Reference stress subfield." + + referenceStrain = pythia.pyre.inventory.facility("reference_strain", family="auxiliary_subfield", factory=Subfield) + referenceStrain.meta['tip'] = "Reference strain subfield." + # PUBLIC METHODS ///////////////////////////////////////////////////// def __init__(self, name="auxfieldsisotropiclinearporoelasticity"): diff --git a/pylith/problems/SolnDispLagrange.py b/pylith/problems/SolnDispLagrange.py index 6458a197bd..013cd3b1cf 100644 --- a/pylith/problems/SolnDispLagrange.py +++ b/pylith/problems/SolnDispLagrange.py @@ -40,7 +40,7 @@ class SolnDispLagrange(PetscComponent): displacement.meta['tip'] = "Displacement subfield." from .SubfieldLagrangeFault import SubfieldLagrangeFault - lagrangeFault = pythia.pyre.inventory.facility("lagrange_fault", family="soln_subfield", factory=SubfieldLagrangeFault) + lagrangeFault = pythia.pyre.inventory.facility("lagrange_multiplier_fault", family="soln_subfield", factory=SubfieldLagrangeFault) lagrangeFault.meta['tip'] = "Fault Lagrange multiplier subfield." # PUBLIC METHODS ///////////////////////////////////////////////////// @@ -57,7 +57,7 @@ def _configure(self): def components(self): """Order of facilities in Inventory is ambiguous, so overwrite - components() to insure order is [displacement, lagrange_fault]. + components() to insure order is [displacement, lagrange_multiplier_fault]. """ return [self.displacement, self.lagrangeFault] diff --git a/pylith/problems/SolnDispPresLagrange.py b/pylith/problems/SolnDispPresLagrange.py index d77a983c93..35e0e31b48 100644 --- a/pylith/problems/SolnDispPresLagrange.py +++ b/pylith/problems/SolnDispPresLagrange.py @@ -39,7 +39,7 @@ class SolnDispPresLagrange(PetscComponent): pressure.meta['tip'] = "Pressure subfield." from .SubfieldLagrangeFault import SubfieldLagrangeFault - lagrangeFault = pythia.pyre.inventory.facility("lagrange_fault", family="soln_subfield", factory=SubfieldLagrangeFault) + lagrangeFault = pythia.pyre.inventory.facility("lagrange_multiplier_fault", family="soln_subfield", factory=SubfieldLagrangeFault) lagrangeFault.meta['tip'] = "Fault Lagrange multiplier subfield." def __init__(self, name="solndisppres"): @@ -52,7 +52,7 @@ def _configure(self): def components(self): """Order of facilities in Inventory is ambiguous, so overwrite - components() to insure order is [displacement, pressure, lagrange_fault]. + components() to insure order is [displacement, pressure, lagrange_multiplier_fault]. """ return [self.displacement, self.pressure, self.lagrangeFault] diff --git a/pylith/problems/SolnDispPresTracStrainLagrange.py b/pylith/problems/SolnDispPresTracStrainLagrange.py new file mode 100644 index 0000000000..742c3dc98b --- /dev/null +++ b/pylith/problems/SolnDispPresTracStrainLagrange.py @@ -0,0 +1,84 @@ +# ---------------------------------------------------------------------- +# +# Brad T. Aagaard, U.S. Geological Survey +# Charles A. Williams, GNS Science +# Matthew G. Knepley, University at Buffalo +# +# This code was developed as part of the Computational Infrastructure +# for Geodynamics (http://geodynamics.org). +# +# Copyright (c) 2010-2022 University of California, Davis +# +# See LICENSE.md for license information. +# +# ---------------------------------------------------------------------- + +from pylith.utils.PetscComponent import PetscComponent +from .Solution import Solution as SolutionBase + + +class SolnDispPresTracStrainLagrange(PetscComponent): + """ + Container for solution subfields with displacement, pore pressure, trace strain subfields, and fault Lagrange multiplier subfields. + """ + DOC_CONFIG = { + "cfg": """ + [pylithapp.problem] + solution = pylith.problems.SolnDispPresTracStrainLagrange + """ + } + + import pythia.pyre.inventory + + from .SubfieldDisplacement import SubfieldDisplacement + displacement = pythia.pyre.inventory.facility("displacement", family="soln_subfield", factory=SubfieldDisplacement) + displacement.meta['tip'] = "Displacement subfield." + + from .SubfieldPressure import SubfieldPressure + pressure = pythia.pyre.inventory.facility("pressure", family="soln_subfield", factory=SubfieldPressure) + pressure.meta['tip'] = "Pressure subfield." + + from .SubfieldTraceStrain import SubfieldTraceStrain + traceStrain = pythia.pyre.inventory.facility("trace_strain", family="soln_subfield", factory=SubfieldTraceStrain) + traceStrain.meta['tip'] = "Trace strain subfield." + + from .SubfieldLagrangeFault import SubfieldLagrangeFault + lagrangeFault = pythia.pyre.inventory.facility("lagrange_multiplier_fault", family="soln_subfield", factory=SubfieldLagrangeFault) + lagrangeFault.meta['tip'] = "Fault Lagrange multiplier subfield." + + def __init__(self, name="SolnDispPresTracStrainLagrange"): + """Constructor. + """ + PetscComponent.__init__(self, name, facility="soln_subfields") + + def _configure(self): + PetscComponent._configure(self) + + def components(self): + """Order of facilities in Inventory is ambiguous, so overwrite + components() to insure order is [displacement, pressure, trace_strain, lagrange_multiplier_fault]. + + """ + return [self.displacement, self.pressure, self.traceStrain, self.lagrangeFault] + + +class Solution(SolutionBase): + """Python solution field with displacement, pressure, trace strain,and Lagrange multiplier subfields. + """ + + import pythia.pyre.inventory + + from .SolutionSubfield import subfieldFactory + subfields = pythia.pyre.inventory.facilityArray( + "subfields", family="soln_subfields", itemFactory=subfieldFactory, factory=SolnDispPresTracStrainLagrange) + subfields.meta['tip'] = "Subfields in solution." + + +# FACTORIES //////////////////////////////////////////////////////////// +def solution(): + """Factory associated with Solution. + """ + return Solution() + + +# End of file diff --git a/pylith/problems/SolnDispPresTracStrainVelPdotTdotLagrange.py b/pylith/problems/SolnDispPresTracStrainVelPdotTdotLagrange.py new file mode 100644 index 0000000000..1762ea43b0 --- /dev/null +++ b/pylith/problems/SolnDispPresTracStrainVelPdotTdotLagrange.py @@ -0,0 +1,107 @@ +# ---------------------------------------------------------------------- +# +# Brad T. Aagaard, U.S. Geological Survey +# Charles A. Williams, GNS Science +# Matthew G. Knepley, University at Buffalo +# +# This code was developed as part of the Computational Infrastructure +# for Geodynamics (http://geodynamics.org). +# +# Copyright (c) 2010-2022 University of California, Davis +# +# See LICENSE.md for license information. +# +# ---------------------------------------------------------------------- +# + +# @file pylith/problems/SolnDispPresTracStrainVelPdotTdotLagrange.py +# +# @brief Python subfields container with displacement, pore pressure, and trace strain subfields, along with their time derivatives. + +from pylith.utils.PetscComponent import PetscComponent +from .Solution import Solution as SolutionBase + + +class SolnDispPresTracStrainVelPdotTdotLagrange(PetscComponent): + """ + Python solution field with displacement, pressure, and trace strain subfields, along with their time derivatives, and a lagrange fault. + """ + DOC_CONFIG = { + "cfg": """ + [pylithapp.problem] + solution = pylith.problems.SolnDispPresTracStrainVelPdotTdotLagrange + """ + } + + import pythia.pyre.inventory + + from .SubfieldDisplacement import SubfieldDisplacement + displacement = pythia.pyre.inventory.facility( + "displacement", family="soln_subfield", factory=SubfieldDisplacement) + displacement.meta['tip'] = "Displacement subfield." + + from .SubfieldPressure import SubfieldPressure + pressure = pythia.pyre.inventory.facility( + "pressure", family="soln_subfield", factory=SubfieldPressure) + pressure.meta['tip'] = "Pressure subfield." + + from .SubfieldTraceStrain import SubfieldTraceStrain + traceStrain = pythia.pyre.inventory.facility( + "trace_strain", family="soln_subfield", factory=SubfieldTraceStrain) + traceStrain.meta['tip'] = "Trace strain subfield." + + from .SubfieldVelocity import SubfieldVelocity + velocity = pythia.pyre.inventory.facility( + "velocity", family="soln_subfield", factory=SubfieldVelocity) + velocity.meta['tip'] = "Velocity subfield." + + from .SubfieldPressureDot import SubfieldPressureDot + pressureT = pythia.pyre.inventory.facility( + "pressure_t", family="soln_subfield", factory=SubfieldPressureDot) + pressureT.meta['tip'] = "PressureT subfield." + + from .SubfieldTraceStrainDot import SubfieldTraceStrainDot + traceStrainT = pythia.pyre.inventory.facility( + "trace_strain_t", family="soln_subfield", factory=SubfieldTraceStrainDot) + traceStrainT.meta['tip'] = "TraceStrainT subfield." + + from .SubfieldLagrangeFault import SubfieldLagrangeFault + lagrangeFault = pythia.pyre.inventory.facility("lagrange_multiplier_fault", family="soln_subfield", factory=SubfieldLagrangeFault) + lagrangeFault.meta['tip'] = "Fault Lagrange multiplier subfield." + + def __init__(self, name="solndispprestracstrainvelpdottdotlagrange"): + """Constructor. + """ + PetscComponent.__init__(self, name, facility="soln_subfields") + + def _configure(self): + PetscComponent._configure(self) + + def components(self): + """Order of facilities in Inventory is ambiguous, so overwrite + components() to insure order is [displacement, pressure, trace_strain, velocity, pressure_t, trace_strain_t, lagrange_multiplier_fault]. + + """ + return [self.displacement, self.pressure, self.traceStrain, self.velocity, self.pressureT, self.traceStrainT, self.lagrangeFault] + + +class Solution(SolutionBase): + """Python solution field with displacement, pressure, and trace strain subfields, along with their time derivatives, and a lagrange fault. + """ + + import pythia.pyre.inventory + + from .SolutionSubfield import subfieldFactory + subfields = pythia.pyre.inventory.facilityArray( + "subfields", family="soln_subfields", itemFactory=subfieldFactory, factory=SolnDispPresTracStrainVelPdotTdotLagrange) + subfields.meta['tip'] = "Subfields in solution." + + +# FACTORIES //////////////////////////////////////////////////////////// +def solution(): + """Factory associated with Solution. + """ + return Solution() + + +# End of file \ No newline at end of file diff --git a/pylith/problems/SolnDispVelLagrange.py b/pylith/problems/SolnDispVelLagrange.py index 6c7e6bf2b4..8221fcfde6 100644 --- a/pylith/problems/SolnDispVelLagrange.py +++ b/pylith/problems/SolnDispVelLagrange.py @@ -39,7 +39,7 @@ class SolnDispVelLagrange(PetscComponent): velocity.meta['tip'] = "Velocity subfield." from .SubfieldLagrangeFault import SubfieldLagrangeFault - lagrangeFault = pythia.pyre.inventory.facility("lagrange_fault", family="soln_subfield", factory=SubfieldLagrangeFault) + lagrangeFault = pythia.pyre.inventory.facility("lagrange_multiplier_fault", family="soln_subfield", factory=SubfieldLagrangeFault) lagrangeFault.meta['tip'] = "Fault Lagrange multiplier subfield." def __init__(self, name="solndispvel"): @@ -52,7 +52,7 @@ def _configure(self): def components(self): """Order of facilities in Inventory is ambiguous, so overwrite - components() to insure order is [displacement, velocity, lagrange_fault]. + components() to insure order is [displacement, velocity, lagrange_multiplier_fault]. """ return [self.displacement, self.velocity, self.lagrangeFault] diff --git a/pylith/problems/SubfieldLagrangeFault.py b/pylith/problems/SubfieldLagrangeFault.py index 19ab2e79b7..0b7a5af3ec 100644 --- a/pylith/problems/SubfieldLagrangeFault.py +++ b/pylith/problems/SubfieldLagrangeFault.py @@ -24,7 +24,7 @@ class SubfieldLagrangeFault(SolutionSubfield): """ DOC_CONFIG = { "cfg": """ - [pylithapp.problems.solution.subfields.lagrange_fault] + [pylithapp.problems.solution.subfields.lagrange_multiplier_fault] alias = lagrange_multiplier_fault basis_order = 1 """ diff --git a/pylith/problems/__init__.py b/pylith/problems/__init__.py index 5e20d99796..0ccfbfb6a1 100644 --- a/pylith/problems/__init__.py +++ b/pylith/problems/__init__.py @@ -33,7 +33,9 @@ "SolnDispVelLagrange", "SolnDispVel", "SolnDispPresTracStrain", - "SolnDispPresTracStrainVelPdotTdot", + "SolnDispPresTracStrainLagrange", + "SolnDispPresTracStrainVelPdotTdot", + "SolnDispPresTracStrainVelPdotTdotLagrange", "SolnDispPresVel", "SolutionSubfield", "SubfieldDisplacement", diff --git a/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg b/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg index 40cd56965f..17d78dc9df 100644 --- a/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg +++ b/tests/fullscale/linearelasticity/faults-3d-buried/pylithapp.cfg @@ -65,7 +65,7 @@ solution = pylith.problems.SolnDispLagrange [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 [pylithapp.problem] solution_observers = [domain, bc_ypos] diff --git a/tests/fullscale/linearelasticity/faults-3d/pylithapp.cfg b/tests/fullscale/linearelasticity/faults-3d/pylithapp.cfg index e5a5cfb715..454ac62a3c 100644 --- a/tests/fullscale/linearelasticity/faults-3d/pylithapp.cfg +++ b/tests/fullscale/linearelasticity/faults-3d/pylithapp.cfg @@ -65,7 +65,7 @@ solution = pylith.problems.SolnDispLagrange [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 [pylithapp.problem] solution_observers = [domain, bc_ypos, points] diff --git a/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b0.cfg b/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b0.cfg index 33cd5d4eb2..5c4eb996b0 100644 --- a/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b0.cfg +++ b/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b0.cfg @@ -16,7 +16,7 @@ defaults.quadrature_order = 1 [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 # ---------------------------------------------------------------------- # fault diff --git a/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b1.cfg b/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b1.cfg index 4c551372fb..a9a3f405c9 100644 --- a/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b1.cfg +++ b/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b1.cfg @@ -16,7 +16,7 @@ defaults.quadrature_order = 1 [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 # ---------------------------------------------------------------------- # fault diff --git a/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b2.cfg b/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b2.cfg index 579765f509..3265b9dd06 100644 --- a/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b2.cfg +++ b/tests/fullscale/linearelasticity/greensfns-2d/leftlateral_b2.cfg @@ -16,7 +16,7 @@ defaults.quadrature_order = 2 [pylithapp.problem.solution.subfields] displacement.basis_order = 2 -lagrange_fault.basis_order = 2 +lagrange_multiplier_fault.basis_order = 2 # ---------------------------------------------------------------------- # fault diff --git a/tests/fullscale/linearelasticity/greensfns-2d/opening.cfg b/tests/fullscale/linearelasticity/greensfns-2d/opening.cfg index b5c668709e..43e4a2774e 100644 --- a/tests/fullscale/linearelasticity/greensfns-2d/opening.cfg +++ b/tests/fullscale/linearelasticity/greensfns-2d/opening.cfg @@ -16,7 +16,7 @@ defaults.quadrature_order = 1 [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 # ---------------------------------------------------------------------- # fault diff --git a/tests/fullscale/linearelasticity/greensfns-2d/slipthreshold.cfg b/tests/fullscale/linearelasticity/greensfns-2d/slipthreshold.cfg index 0fef2e06bf..40be05aae1 100644 --- a/tests/fullscale/linearelasticity/greensfns-2d/slipthreshold.cfg +++ b/tests/fullscale/linearelasticity/greensfns-2d/slipthreshold.cfg @@ -17,7 +17,7 @@ defaults.quadrature_order = 1 [pylithapp.problem.solution.subfields] displacement.basis_order = 1 -lagrange_fault.basis_order = 1 +lagrange_multiplier_fault.basis_order = 1 # ---------------------------------------------------------------------- # fault diff --git a/tests/manual/2d/faults_2d_buried/faulttest.cfg b/tests/manual/2d/faults_2d_buried/faulttest.cfg index d82b8e5de5..53f3ff205f 100644 --- a/tests/manual/2d/faults_2d_buried/faulttest.cfg +++ b/tests/manual/2d/faults_2d_buried/faulttest.cfg @@ -16,8 +16,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 # ---------------------------------------------------------------------- # fault