Add ability to check scale of material properties when getting values from spatial databases

CHANGES.md

@@ -10,6 +10,7 @@ The version numbers are in the form `MAJOR.MINOR.PATCH`, where major releases in
 * **Added**
   * Default filenames for progress monitor and parameters file are set from the simulation name like
   the other output files.
+  * Consistency check for material properties and scales used in nondimensionalization (currently just the shear modulus).
   * Add section to User Guide on troubleshooting solver issues.
 * **Changed**
   * Switch CI from Azure Pipelines to GitHub Actions.

libsrc/pylith/materials/AuxiliaryFactoryElastic.cc

@@ -55,6 +55,7 @@ pylith::materials::AuxiliaryFactoryElastic::addShearModulus(void) {
     description.componentNames[0] = subfieldName;
     description.scale = pressureScale;
     description.validator = pylith::topology::FieldQuery::validatorNonnegative;
+    description.validatorTolerance = 25.0;
 
     _field->subfieldAdd(description, getSubfieldDiscretization(subfieldName));
     pylith::materials::Query::shearModulusFromVM(subfieldName, this);

libsrc/pylith/materials/AuxiliaryFactoryElasticity.cc

@@ -72,7 +72,11 @@ pylith::materials::AuxiliaryFactoryElasticity::addBodyForce(void) {
     PYLITH_JOURNAL_DEBUG("addBodyForce(void)");
 
     const char* subfieldName = "body_force";
-    const char* componentNames[3] = { "body_force_x", "body_force_y", "body_force_z" };
+    const char* componentNames[3] = {
+        "body_force_x",
+        "body_force_y",
+        "body_force_z",
+    };
 
     const PylithReal forceScale = _normalizer->getPressureScale() / _normalizer->getLengthScale();
 

libsrc/pylith/materials/AuxiliaryFactoryPoroelastic.cc

@@ -39,7 +39,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::~AuxiliaryFactoryPoroelastic(voi
 // ----------------------------------------------------------------------------
 // Add isotropic permeability subfield to auxiliary fields.
 void
-pylith::materials::AuxiliaryFactoryPoroelastic::addIsotropicPermeability(void) { // isotropicPermeablity
+pylith::materials::AuxiliaryFactoryPoroelastic::addIsotropicPermeability(void) {
     PYLITH_METHOD_BEGIN;
     PYLITH_JOURNAL_DEBUG("addIsotropicPermeability(void)");
 
@@ -56,7 +56,8 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addIsotropicPermeability(void) {
     description.componentNames.resize(1);
     description.componentNames[0] = subfieldName;
     description.scale = permeabilityScale;
-    description.validator = NULL;
+    // description.validator = pylith::topology::FieldQuery::validatorScale;
+    // description.validatorTolerance = 10.0;
 
     _field->subfieldAdd(description, getSubfieldDiscretization(subfieldName));
     this->setSubfieldQuery(subfieldName);
@@ -68,7 +69,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addIsotropicPermeability(void) {
 // ----------------------------------------------------------------------------
 // Add isotropic permeability subfield to auxiliary fields.
 void
-pylith::materials::AuxiliaryFactoryPoroelastic::addTensorPermeability(void) { // tensorPermeablity
+pylith::materials::AuxiliaryFactoryPoroelastic::addTensorPermeability(void) {
     PYLITH_METHOD_BEGIN;
     PYLITH_JOURNAL_DEBUG("addTensorPermeability(void)");
 
@@ -95,7 +96,8 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addTensorPermeability(void) { //
         description.componentNames[i] = componentNames[i];
     } // for
     description.scale = permeabilityScale;
-    description.validator = NULL;
+    // description.validator = pylith::topology::FieldQuery::validatorScale;
+    // description.validatorTolerance = 10.0;
 
     _field->subfieldAdd(description, getSubfieldDiscretization(subfieldName));
     this->setSubfieldQuery(subfieldName);
@@ -107,7 +109,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addTensorPermeability(void) { //
 // --------------------------------------------------------------------
 // Add drained bulk modulus subfield to auxiliary fields.
 void
-pylith::materials::AuxiliaryFactoryPoroelastic::addDrainedBulkModulus(void) { // DrainedBulkModulus
+pylith::materials::AuxiliaryFactoryPoroelastic::addDrainedBulkModulus(void) {
     PYLITH_METHOD_BEGIN;
     PYLITH_JOURNAL_DEBUG("addDrainedBulkModulus(void)");
 
@@ -123,6 +125,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addDrainedBulkModulus(void) { //
     description.componentNames[0] = subfieldName;
     description.scale = pressureScale;
     description.validator = pylith::topology::FieldQuery::validatorPositive;
+    description.validatorTolerance = 25.0;
 
     _field->subfieldAdd(description, getSubfieldDiscretization(subfieldName));
     this->setSubfieldQuery(subfieldName);
@@ -134,7 +137,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addDrainedBulkModulus(void) { //
 // ---------------------------------------------------------------------
 // Add biot coefficient subfield to auxiliary fields.
 void
-pylith::materials::AuxiliaryFactoryPoroelastic::addBiotCoefficient(void) { // biotCoefficient
+pylith::materials::AuxiliaryFactoryPoroelastic::addBiotCoefficient(void) {
     PYLITH_METHOD_BEGIN;
     PYLITH_JOURNAL_DEBUG("addBiotCoefficient(void)");
 
@@ -149,6 +152,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addBiotCoefficient(void) { // bi
     description.componentNames[0] = subfieldName;
     description.scale = 1.0;
     description.validator = pylith::topology::FieldQuery::validatorPositive;
+    description.validatorTolerance = 10.0;
 
     _field->subfieldAdd(description, getSubfieldDiscretization(subfieldName));
     this->setSubfieldQuery(subfieldName);
@@ -160,7 +164,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addBiotCoefficient(void) { // bi
 // ---------------------------------------------------------------------
 // Add biot modulus subfield to auxiliary fields.
 void
-pylith::materials::AuxiliaryFactoryPoroelastic::addBiotModulus(void) { // biotCoefficient
+pylith::materials::AuxiliaryFactoryPoroelastic::addBiotModulus(void) {
     PYLITH_METHOD_BEGIN;
     PYLITH_JOURNAL_DEBUG("addBiotModulus(void)");
 
@@ -176,6 +180,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addBiotModulus(void) { // biotCo
     description.componentNames[0] = subfieldName;
     description.scale = pressureScale;
     description.validator = pylith::topology::FieldQuery::validatorPositive;
+    description.validatorTolerance = 25.0;
 
     _field->subfieldAdd(description, getSubfieldDiscretization(subfieldName));
     pylith::materials::Query::biotModulusFromInput(subfieldName, this);
@@ -187,7 +192,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addBiotModulus(void) { // biotCo
 // ----------------------------------------------------------------------
 // Add reference stress subfield to auxiliary fields.
 void
-pylith::materials::AuxiliaryFactoryPoroelastic::addReferenceStress(void) { // referenceStress
+pylith::materials::AuxiliaryFactoryPoroelastic::addReferenceStress(void) {
     PYLITH_METHOD_BEGIN;
     PYLITH_JOURNAL_DEBUG("addReferenceStress(void)");
 
@@ -225,7 +230,7 @@ pylith::materials::AuxiliaryFactoryPoroelastic::addReferenceStress(void) { // re
 // ----------------------------------------------------------------------
 // Add reference strain subfield to auxiliary fields.
 void
-pylith::materials::AuxiliaryFactoryPoroelastic::addReferenceStrain(void) { // addReferenceStrain
+pylith::materials::AuxiliaryFactoryPoroelastic::addReferenceStrain(void) {
     PYLITH_METHOD_BEGIN;
     PYLITH_JOURNAL_DEBUG("addRefrenceStrain(void)");
 

libsrc/pylith/materials/AuxiliaryFactoryPoroelasticity.cc

@@ -116,7 +116,7 @@ pylith::materials::AuxiliaryFactoryPoroelasticity::addPorosity(void) { // porosi
     PYLITH_JOURNAL_DEBUG("addPorosity(void)");
 
     const char* subfieldName = "porosity";
-    const PylithReal noScale = 1;
+    const PylithReal noScale = 1.0;
 
     pylith::topology::Field::Description description;
     description.label = subfieldName;

libsrc/pylith/materials/AuxiliaryFactoryPoroelasticity.hh

@@ -10,11 +10,11 @@
 #pragma once
 
 #include "pylith/materials/materialsfwd.hh" // forward declarations
-#include "pylith/feassemble/AuxiliaryFactory.hh" // ISA AuxiliaryFactory
+#include "pylith/materials/AuxiliaryFactoryElasticity.hh" // ISA AuxiliaryFactory
 
 #include "spatialdata/spatialdb/spatialdbfwd.hh" // USES GravityField
 
-class pylith::materials::AuxiliaryFactoryPoroelasticity : public pylith::feassemble::AuxiliaryFactory {
+class pylith::materials::AuxiliaryFactoryPoroelasticity : public pylith::materials::AuxiliaryFactoryElasticity {
     friend class TestAuxiliaryFactoryPoroelasticity; // unit testing
 
     // PUBLIC METHODS //////////////////////////////////////////////////////////////////////////////////////////////////

libsrc/pylith/materials/AuxiliaryFactoryViscoelastic.cc

@@ -114,8 +114,11 @@ pylith::materials::AuxiliaryFactoryViscoelastic::addShearModulusRatioGeneralized
     description.vectorFieldType = pylith::topology::Field::OTHER;
     description.numComponents = 3;
     description.componentNames.resize(3);
-    const char* componentNames[3] = { "shear_modulus_ratio_1", "shear_modulus_ratio_2",
-                                      "shear_modulus_ratio_3" };
+    const char* componentNames[3] = {
+        "shear_modulus_ratio_1",
+        "shear_modulus_ratio_2",
+        "shear_modulus_ratio_3",
+    };
     for (int i = 0; i < 3; ++i) {
         description.componentNames[i] = componentNames[i];
     } // for

libsrc/pylith/materials/Query.cc

@@ -222,9 +222,9 @@ pylith::materials::_Query::vmToShearModulus(PylithScalar valueSubfield[],
 
     std::ostringstream msg;
     if (density <= 0) {
-        msg << "Found negative density (" << density << ").";
+        msg << "Found non-positive density (" << density << ").";
     } // if
-    if (vs <= 0) {
+    if (vs < 0) {
         msg << "Found negative shear wave speed (" << vs << ").";
     } // if
 
@@ -257,13 +257,13 @@ pylith::materials::_Query::vmToBulkModulus(PylithScalar valueSubfield[],
 
     std::ostringstream msg;
     if (density <= 0) {
-        msg << "Found nonpositive density (" << density << ").";
+        msg << "Found non-positive density (" << density << ").";
     } // if
     if (vs < 0) {
         msg << "Found negative shear wave speed (" << vs << ").";
     } // if
     if (vp <= 0) {
-        msg << "Found nonpositive dilatational wave speed (" << vp << ").";
+        msg << "Found non-positive dilatational wave speed (" << vp << ").";
     } // if
 
     PYLITH_METHOD_RETURN(msg.str());
@@ -296,13 +296,13 @@ pylith::materials::_Query::vmToMaxwellTime(PylithScalar valueSubfield[],
 
     std::ostringstream msg;
     if (density <= 0) {
-        msg << "Found negative density (" << density << ").";
+        msg << "Found non-positive density (" << density << ").";
     } // if
     if (vs <= 0) {
-        msg << "Found negative shear wave speed (" << vs << ").";
+        msg << "Found non-positive shear wave speed (" << vs << ").";
     } // if
     if (viscosity <= 0) {
-        msg << "Found nonpositive viscosity (" << viscosity << ").";
+        msg << "Found non-positive viscosity (" << viscosity << ").";
     } // if
 
     PYLITH_METHOD_RETURN(msg.str());
@@ -351,28 +351,28 @@ pylith::materials::_Query::vmToGeneralizedMaxwellTimes(PylithScalar valueSubfiel
 
     std::ostringstream msg;
     if (density <= 0) {
-        msg << "Found negative density (" << density << ").";
+        msg << "Found non-positive density (" << density << ").";
     } // if
     if (vs <= 0) {
-        msg << "Found negative shear wave speed (" << vs << ").";
+        msg << "Found non-positive shear wave speed (" << vs << ").";
     } // if
     if (viscosity1 <= 0) {
-        msg << "Found nonpositive viscosity 1 (" << viscosity1 << ").";
+        msg << "Found non-positive viscosity 1 (" << viscosity1 << ").";
     } // if
     if (viscosity2 <= 0) {
-        msg << "Found nonpositive viscosity 2 (" << viscosity2 << ").";
+        msg << "Found non-positive viscosity 2 (" << viscosity2 << ").";
     } // if
     if (viscosity3 <= 0) {
-        msg << "Found nonpositive viscosity 3 (" << viscosity3 << ").";
+        msg << "Found non-positive viscosity 3 (" << viscosity3 << ").";
     } // if
     if (shearModulusRatio1 <= 0) {
-        msg << "Found negative shear modulus ratio 1 (" << shearModulusRatio1 << ").";
+        msg << "Found non-positive shear modulus ratio 1 (" << shearModulusRatio1 << ").";
     } // if
     if (shearModulusRatio2 <= 0) {
-        msg << "Found negative shear modulus ratio 2 (" << shearModulusRatio2 << ").";
+        msg << "Found non-positive shear modulus ratio 2 (" << shearModulusRatio2 << ").";
     } // if
     if (shearModulusRatio3 <= 0) {
-        msg << "Found negative shear modulus ratio 3 (" << shearModulusRatio3 << ").";
+        msg << "Found non-positive shear modulus ratio 3 (" << shearModulusRatio3 << ").";
     } // if
 
     const double ratioSum = shearModulusRatio1 + shearModulusRatio2 + shearModulusRatio3;
@@ -409,13 +409,13 @@ pylith::materials::_Query::vmToGeneralizedMaxwellShearModulusRatios(PylithScalar
 
     std::ostringstream msg;
     if (shearModulusRatio1 <= 0) {
-        msg << "Found negative shear modulus ratio 1 (" << shearModulusRatio1 << ").";
+        msg << "Found non-positive shear modulus ratio 1 (" << shearModulusRatio1 << ").";
     } // if
     if (shearModulusRatio2 <= 0) {
-        msg << "Found negative shear modulus ratio 2 (" << shearModulusRatio2 << ").";
+        msg << "Found non-positive shear modulus ratio 2 (" << shearModulusRatio2 << ").";
     } // if
     if (shearModulusRatio3 <= 0) {
-        msg << "Found negative shear modulus ratio 3 (" << shearModulusRatio3 << ").";
+        msg << "Found non-positive shear modulus ratio 3 (" << shearModulusRatio3 << ").";
     } // if
     const double ratioSum = shearModulusRatio1 + shearModulusRatio2 + shearModulusRatio3;
     if (ratioSum > 1) {
@@ -498,7 +498,7 @@ pylith::materials::_Query::inputToBiotModulus(PylithScalar valueSubfield[],
         msg << "Found negative porosity (" << porosity << ").";
     } // if
     if (biot_coefficient <= 0) {
-        msg << "Found negative biot coefficient (" << biot_coefficient << ").";
+        msg << "Found non-positive biot coefficient (" << biot_coefficient << ").";
     } // if
 
     // Debug

libsrc/pylith/topology/FieldBase.hh

@@ -46,7 +46,9 @@ public:
 public:
 
     /// Function prototype for validator functions.
-    typedef const char* (*validatorfn_type)(const PylithReal);
+    typedef std::string (*validatorfn_type)(const PylithReal /*value */,
+                                            const PylithReal /* scale */,
+                                            const PylithReal /* tolerance */);
 
     // PUBLIC STRUCTS ///////////////////////////////////////////////////////
 public:
@@ -58,6 +60,7 @@ public:
         pylith::string_vector componentNames; ///< Names of components.
         size_t numComponents; ///< Number of components.
         PylithReal scale; ///< Dimension scale associated with values.
+        PylithReal validatorTolerance; ///< Tolerance relative to scale for validation.
         validatorfn_type validator; ///< Validator for values in field;
         bool hasHistory; ///< Has subfields with history, i.e., state variables.
         size_t historySize; ///< Number of points in time history (currently only).
@@ -68,6 +71,7 @@ public:
             vectorFieldType(OTHER),
             numComponents(0),
             scale(1.0),
+            validatorTolerance(0.0),
             validator(NULL),
             hasHistory(false),
             historySize(0) {}
@@ -79,6 +83,7 @@ public:
                     const size_t numComponentsValue=0,
                     const VectorFieldEnum vectorFieldTypeValue=SCALAR,
                     const PylithReal scaleValue=1.0,
+                    const PylithReal validatorToleranceValue=0.0,
                     const validatorfn_type validatorValue=NULL,
                     bool isFaultOnlyValue=false,
                     bool hasHistoryValue=false,
@@ -89,6 +94,7 @@ public:
             componentNames(componentNamesValue),
             numComponents(numComponentsValue),
             scale(scaleValue),
+            validatorTolerance(validatorToleranceValue),
             validator(validatorValue),
             hasHistory(hasHistoryValue),
             historySize(historySizeValue) {}