Merge branch 'PassVariablesPrevInMplCalls' into 'master'

Pass variables prev to MPL Property::value() calls.

See merge request ogs/ogs!3102

MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp

@@ -34,10 +34,21 @@ void PorosityFromMassBalance::checkScale() const
     }
 }
 
+PropertyDataType PorosityFromMassBalance::value(
+    VariableArray const& /*variable_array*/,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    OGS_FATAL(
+        "PorosityFromMassBalance value call requires previous time step "
+        "values.");
+}
+
 PropertyDataType PorosityFromMassBalance::value(
     VariableArray const& variable_array,
-    ParameterLib::SpatialPosition const& pos,
-    double const t, double const dt) const
+    VariableArray const& variable_array_prev,
+    ParameterLib::SpatialPosition const& pos, double const t,
+    double const dt) const
 {
     double const beta_SR = std::get<double>(
         variable_array[static_cast<int>(Variable::grain_compressibility)]);
@@ -46,17 +57,24 @@ PropertyDataType PorosityFromMassBalance::value(
             ->property(PropertyType::biot_coefficient)
             .template value<double>(variable_array, pos, t, dt);
 
-    double const e_dot = std::get<double>(
-        variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
+    double const e = std::get<double>(
+        variable_array[static_cast<int>(Variable::volumetric_strain)]);
+    double const e_prev = std::get<double>(
+        variable_array_prev[static_cast<int>(Variable::volumetric_strain)]);
+    double const delta_e = e - e_prev;
 
-    double const p_eff_dot = std::get<double>(variable_array[static_cast<int>(
-        Variable::effective_pore_pressure_rate)]);
+    double const p_eff = std::get<double>(
+        variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
+    double const p_eff_prev =
+        std::get<double>(variable_array_prev[static_cast<int>(
+            Variable::effective_pore_pressure)]);
+    double const delta_p_eff = p_eff - p_eff_prev;
 
-    double const phi = std::get<double>(
-        variable_array[static_cast<int>(Variable::porosity)]);
+    double const phi_prev = std::get<double>(
+        variable_array_prev[static_cast<int>(Variable::porosity)]);
 
-    double const w = dt * (e_dot + p_eff_dot * beta_SR);
-    return std::clamp((phi + alpha_b * w) / (1 + w), phi_min_, phi_max_);
+    double const w = delta_e + delta_p_eff * beta_SR;
+    return std::clamp((phi_prev + alpha_b * w) / (1 + w), phi_min_, phi_max_);
 }
 
 PropertyDataType PorosityFromMassBalance::dValue(

MaterialLib/MPL/Properties/PorosityFromMassBalance.h

@@ -51,6 +51,10 @@ public:
     PropertyDataType value(VariableArray const& variable_array,
                            ParameterLib::SpatialPosition const& pos,
                            double const t, double const dt) const override;
+    PropertyDataType value(VariableArray const& variable_array,
+                           VariableArray const& variable_array_prev,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t, double const dt) const override;
     PropertyDataType dValue(VariableArray const& variable_array,
                             Variable const variable,
                             ParameterLib::SpatialPosition const& pos,

MaterialLib/MPL/Properties/SaturationDependentSwelling.cpp

@@ -55,13 +55,14 @@ void SaturationDependentSwelling::checkScale() const
 
 PropertyDataType SaturationDependentSwelling::value(
     VariableArray const& variable_array,
+    VariableArray const& variable_array_prev,
     ParameterLib::SpatialPosition const& pos, double const /*t*/,
     double const dt) const
 {
     auto const S_L = std::get<double>(
         variable_array[static_cast<int>(Variable::liquid_saturation)]);
-    auto const S_L_dot = std::get<double>(
-        variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
+    auto const S_L_prev = std::get<double>(
+        variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
 
     Eigen::Matrix<double, 3, 3> const e =
         local_coordinate_system_ == nullptr
@@ -76,8 +77,6 @@ PropertyDataType SaturationDependentSwelling::value(
         return delta_sigma_sw;   // still being zero.
     }
 
-    double const S_L_prev = S_L - S_L_dot * dt;
-
     double const S_eff = std::clamp((S_L - S_min_) / (S_max_ - S_min_), 0., 1.);
     double const S_eff_prev =
         std::clamp((S_L_prev - S_min_) / (S_max_ - S_min_), 0., 1.);
@@ -100,9 +99,10 @@ PropertyDataType SaturationDependentSwelling::value(
 }
 
 PropertyDataType SaturationDependentSwelling::dValue(
-    VariableArray const& variable_array, Variable const primary_variable,
+    VariableArray const& variable_array,
+    VariableArray const& variable_array_prev, Variable const primary_variable,
     ParameterLib::SpatialPosition const& pos, double const /*t*/,
-    double const dt) const
+    double const /*dt*/) const
 {
     (void)primary_variable;
     assert((primary_variable == Variable::liquid_saturation) &&
@@ -111,8 +111,8 @@ PropertyDataType SaturationDependentSwelling::dValue(
 
     auto const S_L = std::get<double>(
         variable_array[static_cast<int>(Variable::liquid_saturation)]);
-    auto const S_L_dot = std::get<double>(
-        variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
+    auto const S_L_prev = std::get<double>(
+        variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
 
     Eigen::Matrix<double, 3, 3> const e =
         local_coordinate_system_ == nullptr
@@ -127,8 +127,6 @@ PropertyDataType SaturationDependentSwelling::dValue(
         return delta_sigma_sw;   // still being zero.
     }
 
-    double const S_L_prev = S_L - S_L_dot * dt;
-
     double const S_eff = std::clamp((S_L - S_min_) / (S_max_ - S_min_), 0., 1.);
     double const S_eff_prev =
         std::clamp((S_L_prev - S_min_) / (S_max_ - S_min_), 0., 1.);

MaterialLib/MPL/Properties/SaturationDependentSwelling.h

@@ -56,9 +56,11 @@ public:
     void checkScale() const override;
 
     PropertyDataType value(VariableArray const& variable_array,
+                           VariableArray const& variable_array_prev,
                            ParameterLib::SpatialPosition const& pos,
                            double const t, double const dt) const override;
     PropertyDataType dValue(VariableArray const& variable_array,
+                            VariableArray const& variable_array_prev,
                             Variable const variable,
                             ParameterLib::SpatialPosition const& pos,
                             double const t, double const dt) const override;

MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.cpp

@@ -24,10 +24,21 @@ LinearSaturationSwellingStress::LinearSaturationSwellingStress(
     name_ = std::move(name);
 }
 
+PropertyDataType LinearSaturationSwellingStress::value(
+    const VariableArray& /*variable_array*/,
+    const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
+    const double /*dt*/) const
+{
+    OGS_FATAL(
+        "LinearSaturationSwellingStress value call requires previous time step "
+        "values.");
+}
+
 PropertyDataType LinearSaturationSwellingStress::value(
     const VariableArray& variable_array,
+    const VariableArray& variable_array_prev,
     const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
-    const double dt) const
+    const double /*dt*/) const
 {
     // Sl <= S_max is guaranteed by the saturation property or
     // the saturation calculation.
@@ -39,12 +50,10 @@ PropertyDataType LinearSaturationSwellingStress::value(
         return 0.0;
     }
 
-    const double dS =
-        dt *
-        std::get<double>(
-            variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
+    const double Sl_prev = std::get<double>(
+        variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
 
-    return coefficient_ * dS;
+    return coefficient_ * (Sl - Sl_prev);
 }
 
 PropertyDataType LinearSaturationSwellingStress::dValue(

MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h

@@ -77,8 +77,14 @@ public:
         }
     }
 
+    PropertyDataType value(VariableArray const& variable_array,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t,
+                           double const dt) const override;
+
     /// \return \f$\Delta  {{\sigma}}^{\text{sw}} \f$.
     PropertyDataType value(VariableArray const& variable_array,
+                           VariableArray const& variable_array_prev,
                            ParameterLib::SpatialPosition const& pos,
                            double const t,
                            double const dt) const override;

MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp

@@ -36,6 +36,7 @@ void TransportPorosityFromMassBalance::checkScale() const
 
 PropertyDataType TransportPorosityFromMassBalance::value(
     VariableArray const& variable_array,
+    VariableArray const& variable_array_prev,
     ParameterLib::SpatialPosition const& pos, double const t,
     double const dt) const
 {
@@ -46,22 +47,39 @@ PropertyDataType TransportPorosityFromMassBalance::value(
             ->property(PropertyType::biot_coefficient)
             .template value<double>(variable_array, pos, t, dt);
 
-    double const e_dot = std::get<double>(
-        variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
+    double const e = std::get<double>(
+        variable_array[static_cast<int>(Variable::volumetric_strain)]);
+    double const e_prev = std::get<double>(
+        variable_array_prev[static_cast<int>(Variable::volumetric_strain)]);
+    double const delta_e = e - e_prev;
 
-    double const p_eff_dot = std::get<double>(variable_array[static_cast<int>(
-        Variable::effective_pore_pressure_rate)]);
+    double const p_eff = std::get<double>(
+        variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
+    double const p_eff_prev =
+        std::get<double>(variable_array_prev[static_cast<int>(
+            Variable::effective_pore_pressure)]);
+    double const delta_p_eff = p_eff - p_eff_prev;
 
     double const phi =
         std::get<double>(variable_array[static_cast<int>(Variable::porosity)]);
 
     double const phi_tr_prev = std::get<double>(
-        variable_array[static_cast<int>(Variable::transport_porosity)]);
+        variable_array_prev[static_cast<int>(Variable::transport_porosity)]);
 
-    double const w = dt * (e_dot + p_eff_dot * beta_SR);
+    double const w = delta_e + delta_p_eff * beta_SR;
     return std::clamp(phi_tr_prev + (alpha_b - phi) * w, phi_min_, phi_max_);
 }
 
+PropertyDataType TransportPorosityFromMassBalance::value(
+    VariableArray const& /*variable_array*/,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    OGS_FATAL(
+        "TransportPorosityFromMassBalance value call requires previous time "
+        "step values.");
+}
+
 PropertyDataType TransportPorosityFromMassBalance::dValue(
     VariableArray const& /*variable_array*/,
     Variable const /*primary_variable*/,

MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.h

@@ -48,6 +48,11 @@ public:
     }
 
     PropertyDataType value(VariableArray const& variable_array,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t,
+                           double const dt) const override;
+    PropertyDataType value(VariableArray const& variable_array,
+                           VariableArray const& variable_array_prev,
                            ParameterLib::SpatialPosition const& pos,
                            double const t, double const dt) const override;
     PropertyDataType dValue(VariableArray const& variable_array,

MaterialLib/MPL/Property.cpp

@@ -74,18 +74,23 @@ PropertyDataType Property::value() const
     return value_;
 }
 
-/// The default implementation of this method only returns the property value
-/// without altering it.
 PropertyDataType Property::value(VariableArray const& /*variable_array*/,
+                                 VariableArray const& /*variable_array_prev*/,
                                  ParameterLib::SpatialPosition const& /*pos*/,
                                  double const /*t*/, double const /*dt*/) const
 {
     return value_;
 }
 
-/// The default implementation of this method only returns the
-/// property value derivative without altering it.
+PropertyDataType Property::value(VariableArray const& variable_array,
+                                 ParameterLib::SpatialPosition const& pos,
+                                 double const t, double const dt) const
+{
+    return value(variable_array, VariableArray{}, pos, t, dt);
+}
+
 PropertyDataType Property::dValue(VariableArray const& /*variable_array*/,
+                                  VariableArray const& /*variable_array_prev*/,
                                   Variable const /*variable*/,
                                   ParameterLib::SpatialPosition const& /*pos*/,
                                   double const /*t*/, double const /*dt*/) const
@@ -93,6 +98,16 @@ PropertyDataType Property::dValue(VariableArray const& /*variable_array*/,
     return dvalue_;
 }
 
+/// The default implementation of this method only returns the
+/// property value derivative without altering it.
+PropertyDataType Property::dValue(VariableArray const& variable_array,
+                                  Variable const variable,
+                                  ParameterLib::SpatialPosition const& pos,
+                                  double const t, double const dt) const
+{
+    return dValue(variable_array, VariableArray{}, variable, pos, t, dt);
+}
+
 /// Default implementation: 2nd derivative of any constant property is zero.
 PropertyDataType Property::d2Value(VariableArray const& /*variable_array*/,
                                    Variable const /*variable*/,

MaterialLib/MPL/Property.h

@@ -54,13 +54,30 @@ public:
     /// This virtual method simply returns the private value_ attribute without
     /// changing it.
     virtual PropertyDataType value() const;
-    /// This virtual method will compute the property value based on the primary
-    /// variables that are passed as arguments.
+    /// This virtual method will compute the property value based on the
+    /// variables that are passed as arguments and the variables from the
+    /// previous time step.
     virtual PropertyDataType value(VariableArray const& variable_array,
+                                   VariableArray const& variable_array_prev,
                                    ParameterLib::SpatialPosition const& pos,
                                    double const t, double const dt) const;
-    /// This virtual method will compute the derivative of a property
-    /// with respect to the given variable pv.
+    /// This virtual method will compute the property value based on the
+    /// variables that are passed as arguments with the default implementation
+    /// using empty variables array for the previous time step.
+    virtual PropertyDataType value(VariableArray const& variable_array,
+                                   ParameterLib::SpatialPosition const& pos,
+                                   double const t, double const dt) const;
+    /// This virtual method will compute the property derivative value based on
+    /// the variables that are passed as arguments and the variables from the
+    /// previous time step.
+    virtual PropertyDataType dValue(VariableArray const& variable_array,
+                                    VariableArray const& variable_array_prev,
+                                    Variable const variable,
+                                    ParameterLib::SpatialPosition const& pos,
+                                    double const t, double const dt) const;
+    /// This virtual method will compute the property derivative value based on
+    /// the variables that are passed as arguments with the default
+    /// implementation using empty variables array for the previous time step.
     virtual PropertyDataType dValue(VariableArray const& variable_array,
                                     Variable const variable,
                                     ParameterLib::SpatialPosition const& pos,
@@ -116,6 +133,29 @@ public:
         }
     }
 
+    template <typename T>
+    T value(VariableArray const& variable_array,
+            VariableArray const& variable_array_prev,
+            ParameterLib::SpatialPosition const& pos, double const t,
+            double const dt) const
+    {
+        try
+        {
+            return std::get<T>(
+                value(variable_array, variable_array_prev, pos, t, dt));
+        }
+        catch (std::bad_variant_access const&)
+        {
+            OGS_FATAL(
+                "The value of {:s} is not of the requested type '{:s}' but a "
+                "{:s}.",
+                description(),
+                typeid(T).name(),
+                property_data_type_names_[value(variable_array,
+                                                variable_array_prev, pos, t, dt)
+                                              .index()]);
+        }
+    }
     template <typename T>
     T value(VariableArray const& variable_array,
             ParameterLib::SpatialPosition const& pos, double const t,
@@ -137,6 +177,28 @@ public:
         }
     }
 
+    template <typename T>
+    T dValue(VariableArray const& variable_array,
+             VariableArray const& variable_array_prev, Variable const variable,
+             ParameterLib::SpatialPosition const& pos, double const t,
+             double const dt) const
+    {
+        try
+        {
+            return std::get<T>(dValue(variable_array, variable_array_prev,
+                                      variable, pos, t, dt));
+        }
+        catch (std::bad_variant_access const&)
+        {
+            OGS_FATAL(
+                "The first derivative value of {:s} is not of the requested "
+                "type '{:s}' but a {:s}.",
+                description(),
+                typeid(T).name(),
+                property_data_type_names_
+                    [dValue(variable_array, variable, pos, t, dt).index()]);
+        }
+    }
     template <typename T>
     T dValue(VariableArray const& variable_array, Variable const variable,
              ParameterLib::SpatialPosition const& pos, double const t,

MaterialLib/MPL/VariableType.h

@@ -46,10 +46,9 @@ enum class Variable : int
     concentration,
     density,
     displacement,
-    effective_pore_pressure_rate,
+    effective_pore_pressure,
     grain_compressibility,
     liquid_saturation,
-    liquid_saturation_rate,
     phase_pressure,
     porosity,
     solid_grain_pressure,
@@ -57,7 +56,7 @@ enum class Variable : int
     stress,
     temperature,
     transport_porosity,
-    volumetric_strain_rate,
+    volumetric_strain,
     number_of_variables
 };
 

Tests/MaterialLib/TestMPLLinearSaturationSwellingStress.cpp

@@ -37,7 +37,7 @@ TEST(MaterialPropertyLib, LinearSaturationSwellingStress)
         "    <reference_saturation> 0.65 </reference_saturation>"
         "</property>";
 
-    auto const swelling_stress_rate =
+    auto const swelling_stress_increment =
         createLinearSaturationSwellingStressModel(xml);
 
     double const coefficient = 1.0e+6;
@@ -49,20 +49,21 @@ TEST(MaterialPropertyLib, LinearSaturationSwellingStress)
     const double dS = S1 - S0;
 
     MaterialPropertyLib::VariableArray variable_array;
-    variable_array[static_cast<int>(
-        MaterialPropertyLib::Variable::liquid_saturation_rate)] = dS / dt;
+    MaterialPropertyLib::VariableArray variable_array_prev;
     variable_array[static_cast<int>(
         MaterialPropertyLib::Variable::liquid_saturation)] = S1;
+    variable_array_prev[static_cast<int>(
+        MaterialPropertyLib::Variable::liquid_saturation)] = S0;
     ParameterLib::SpatialPosition const pos;
     double const time = std::numeric_limits<double>::quiet_NaN();
     double d_sw_expected = coefficient * dS;
-    double d_sw = std::get<double>(
-        swelling_stress_rate->value(variable_array, pos, time, dt));
+    double d_sw = std::get<double>(swelling_stress_increment->value(
+        variable_array, variable_array_prev, pos, time, dt));
     ASSERT_LE(std::fabs(d_sw_expected - d_sw), 1e-19)
-        << "for expected swelling stress rate" << d_sw_expected
-        << " for computed swelling stress rate." << d_sw_expected;
+        << "for expected swelling stress increment" << d_sw_expected
+        << " for computed swelling stress increment." << d_sw_expected;
 
-    double dsw_dS = std::get<double>(swelling_stress_rate->dValue(
+    double dsw_dS = std::get<double>(swelling_stress_increment->dValue(
         variable_array, MaterialPropertyLib::Variable::liquid_saturation, pos,
         time, dt));
     double dsw_dS_expected = coefficient;
@@ -74,15 +75,17 @@ TEST(MaterialPropertyLib, LinearSaturationSwellingStress)
     double const S2 = 0.3;
     variable_array[static_cast<int>(
         MaterialPropertyLib::Variable::liquid_saturation)] = S2;
+    variable_array_prev[static_cast<int>(
+        MaterialPropertyLib::Variable::liquid_saturation)] = S1;
 
-    d_sw = std::get<double>(
-        swelling_stress_rate->value(variable_array, pos, time, dt));
+    d_sw = std::get<double>(swelling_stress_increment->value(
+        variable_array, variable_array_prev, pos, time, dt));
     d_sw_expected = 0.0;
     ASSERT_LE(std::fabs(d_sw_expected - d_sw), 1e-19)
-        << "for expected swelling stress rate" << d_sw_expected
-        << " for computed swelling stress rate." << d_sw_expected;
+        << "for expected swelling stress increment" << d_sw_expected
+        << " for computed swelling stress increment." << d_sw_expected;
 
-    dsw_dS = std::get<double>(swelling_stress_rate->dValue(
+    dsw_dS = std::get<double>(swelling_stress_increment->dValue(
         variable_array, MaterialPropertyLib::Variable::liquid_saturation, pos,
         time, dt));
     dsw_dS_expected = 0.0;