[MPL] Generalise van-Genuchten relation

Making the second (saturation) exponent independent of the first one the
pressure exponent.

Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md

 \copydoc MaterialPropertyLib::SaturationVanGenuchten
+
+Either the `exponent` (which is the pressure exponent \f$m\f$ in the equations)
+must be set and then the saturation exponent is \f$n = 1 / (1 - m)\f$.
+Or both values the `pressure_exponent` and the `saturation_exponent` must be set
+independently.

Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md

-The exponent of the van Genuchten saturation function.
+The "pressure" exponent \f$m\f$ of the van Genuchten saturation function.
+The "saturation" exponent is then set to \f$n = 1 / (1 - m)\f$.
+
+See `pressure_exponent` and `saturation_exponent` to set the values
+independently.

Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_pressure_exponent.md

+The "pressure" exponent \f$m\f$ of the van Genuchten saturation function.
+The "saturation" exponent \f$n\f$ must be set as well.

Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_saturation_exponent.md

+The "saturation" exponent \f$n\f$ of the van Genuchten saturation function.
+The "pressure" exponent \f$m\f$ must be set as well.

MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationVanGenuchten.cpp

@@ -31,14 +31,31 @@ std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten(
     auto const residual_gas_saturation =
         //! \ogs_file_param{properties__property__SaturationVanGenuchten__residual_gas_saturation}
         config.getConfigParameter<double>("residual_gas_saturation");
-    auto const exponent =
-        //! \ogs_file_param{properties__property__SaturationVanGenuchten__exponent}
-        config.getConfigParameter<double>("exponent");
+
+    double pressure_exponent;
+    double saturation_exponent;
+    if (auto const optional_exponent =
+            //! \ogs_file_param{properties__property__SaturationVanGenuchten__exponent}
+        config.getConfigParameterOptional<double>("exponent"))
+    {
+        pressure_exponent = *optional_exponent;
+        saturation_exponent = 1.0 / (1.0 - pressure_exponent);
+    }
+    else
+    {
+        pressure_exponent =
+            //! \ogs_file_param{properties__property__SaturationVanGenuchten__pressure_exponent}
+            config.getConfigParameter<double>("pressure_exponent");
+        saturation_exponent =
+            //! \ogs_file_param{properties__property__SaturationVanGenuchten__saturation_exponent}
+            config.getConfigParameter<double>("saturation_exponent");
+    }
+
     //! \ogs_file_param{properties__property__SaturationVanGenuchten__p_b}
     auto const p_b = config.getConfigParameter<double>("p_b");
 
     return std::make_unique<SaturationVanGenuchten>(
         std::move(property_name), residual_liquid_saturation,
-        residual_gas_saturation, exponent, p_b);
+        residual_gas_saturation, pressure_exponent, saturation_exponent, p_b);
 }
 }  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.cpp

@@ -21,11 +21,13 @@ SaturationVanGenuchten::SaturationVanGenuchten(
     std::string name,
     double const residual_liquid_saturation,
     double const residual_gas_saturation,
-    double const exponent,
+    double const pressure_exponent,
+    double const saturation_exponent,
     double const p_b)
     : S_L_res_(residual_liquid_saturation),
       S_L_max_(1. - residual_gas_saturation),
-      m_(exponent),
+      m_(pressure_exponent),
+      n_(saturation_exponent),
       p_b_(p_b)
 {
     name_ = std::move(name);
@@ -33,10 +35,18 @@ SaturationVanGenuchten::SaturationVanGenuchten(
     if (!(m_ > 0 && m_ < 1))
     {
         OGS_FATAL(
-            "The exponent value m = {:g} of van Genuchten saturation model, is "
-            "out of its range of (0, 1)",
+            "The pressure exponent value m = {:g} of van Genuchten saturation "
+            "model, is out of its range of (0, 1)",
             m_);
     }
+
+    if (n_ < 1)
+    {
+        OGS_FATAL(
+            "The saturation exponent value n = {:g} of van Genuchten "
+            "saturation model, is out of its range of [1, +inf)",
+            n_);
+    }
 }
 
 PropertyDataType SaturationVanGenuchten::value(
@@ -52,8 +62,7 @@ PropertyDataType SaturationVanGenuchten::value(
     }
 
     double const p = p_cap / p_b_;
-    double const n = 1. / (1. - m_);
-    double const p_to_n = std::pow(p, n);
+    double const p_to_n = std::pow(p, n_);
 
     double const S_eff = std::pow(p_to_n + 1., -m_);
     double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
@@ -80,8 +89,7 @@ PropertyDataType SaturationVanGenuchten::dValue(
     }
 
     double const p = p_cap / p_b_;
-    double const n = 1. / (1. - m_);
-    double const p_to_n = std::pow(p, n);
+    double const p_to_n = std::pow(p, n_);
 
     double const S_eff = std::pow(p_to_n + 1., -m_);
     double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
@@ -91,9 +99,8 @@ PropertyDataType SaturationVanGenuchten::dValue(
         return 0.;
     }
 
-    double const dS_eff_dp_cap = -m_ * std::pow(p, n - 1) *
-                                 std::pow(1 + p_to_n, -1. - m_) /
-                                 (p_b_ * (1. - m_));
+    double const dS_eff_dp_cap =
+        -m_ * n_ * p_to_n * S_eff / (p_cap * (p_to_n + 1.));
     return dS_eff_dp_cap * (S_L_max_ - S_L_res_);
 }
 
@@ -117,8 +124,7 @@ PropertyDataType SaturationVanGenuchten::d2Value(
     }
 
     double const p = p_cap / p_b_;
-    double const n = 1. / (1. - m_);
-    double const p_to_n = std::pow(p, n);
+    double const p_to_n = std::pow(p, n_);
 
     double const S_eff = std::pow(p_to_n + 1., -m_);
     double const S = S_eff * S_L_max_ - S_eff * S_L_res_ + S_L_res_;
@@ -129,8 +135,8 @@ PropertyDataType SaturationVanGenuchten::d2Value(
     }
 
     double const d2S_eff_dp_cap2 =
-        m_ * p_to_n * std::pow(p_to_n + 1., -m_ - 2.) * (p_to_n - m_) /
-        (p_cap * p_cap * (m_ - 1.) * (m_ - 1.));
+        m_ * n_ * n_ * p_to_n * S_eff * (p_to_n - m_) /
+        ((p_cap * p_to_n + p_cap) * (p_cap * p_to_n + p_cap));
     return d2S_eff_dp_cap2 * (S_L_max_ - S_L_res_);
 }
 }  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationVanGenuchten.h

@@ -24,13 +24,12 @@ class Component;
  * \f[S_\text{eff}=\frac{S-S_r}{S_{\text{max}}-S_r}.\f]
  * Above, \f$S_r\f$ and \f$S_{\text{max}}\f$ are the residual and the maximum
  * saturations.
- * The exponent \f$m \in (0,1)\f$ and the pressure scaling parameter \f$p_b\f$
- * (it is equal to \f$\rho g/\alpha\f$ in original publication) are given by the
- * user.
- * The scaling parameter \f$p_b\f$ is given in same units as pressure.
+ * The (pressure) exponent \f$m \in (0,1)\f$ and the pressure scaling parameter
+ * \f$p_b\f$ (it is equal to \f$\rho g/\alpha\f$ in original publication) are
+ * given by the user. The scaling parameter \f$p_b\f$ is given in same units as
+ * pressure.
  *
- * In the original work another exponent \f$n\f$ is used, but usually set to
- * \f$n = 1 / (1 - m)\f$, and also in this implementation.
+ * Another (saturaton) exponent \f$n\f$ is usually set to \f$n = 1 / (1 - m)\f$.
  *
  * The saturation is computed from the capillary pressure as follows:
  * \f[S(p_c)=
@@ -49,7 +48,8 @@ public:
     SaturationVanGenuchten(std::string name,
                            double const residual_liquid_saturation,
                            double const residual_gas_saturation,
-                           double const exponent,
+                           double const pressure_exponent,
+                           double const saturation_exponent,
                            double const p_b);
 
     void checkScale() const override
@@ -83,6 +83,7 @@ private:
     double const S_L_res_;
     double const S_L_max_;
     double const m_;
+    double const n_;
     double const p_b_;
 };
 }  // namespace MaterialPropertyLib

Tests/MaterialLib/TestMPLCapillaryPressureVanGenuchten.cpp

@@ -21,21 +21,20 @@ TEST(MaterialPropertyLib, CapillaryPressureVanGenuchten)
 {
     double const residual_liquid_saturation = 0.1;
     double const residual_gas_saturation = 0.05;
-    double const exponent = 0.79;
+    double const pressure_exponent = 0.79;
+    double const saturation_exponent = 1. / (1. - pressure_exponent);
     double const p_b = 10000;
     double const maximum_capillary_pressure = 20000;
 
-    MPL::Property const& pressure =
-        MPL::CapillaryPressureVanGenuchten{"capillary_pressure",
-                                           residual_liquid_saturation,
-                                           residual_gas_saturation,
-                                           exponent,
-                                           p_b,
-                                           maximum_capillary_pressure};
-
-    MPL::Property const& saturation =
-        MPL::SaturationVanGenuchten{"saturation", residual_liquid_saturation,
-                                    residual_gas_saturation, exponent, p_b};
+    MPL::Property const& pressure = MPL::CapillaryPressureVanGenuchten{
+        "capillary_pressure", residual_liquid_saturation,
+        residual_gas_saturation, pressure_exponent,
+        // TODO? saturation_exponent,
+        p_b, maximum_capillary_pressure};
+
+    MPL::Property const& saturation = MPL::SaturationVanGenuchten{
+        "saturation",      residual_liquid_saturation, residual_gas_saturation,
+        pressure_exponent, saturation_exponent,        p_b};
 
     MPL::VariableArray variable_array;
     ParameterLib::SpatialPosition const pos;

Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp

@@ -20,12 +20,13 @@ TEST(MaterialPropertyLib, SaturationVanGenuchten)
 {
     double const residual_liquid_saturation = 0.1;
     double const residual_gas_saturation = 0.05;
-    double const exponent = 0.79;
+    double const pressure_exponent = 0.79;
+    double const saturation_exponent = 1. / (1. - pressure_exponent);
     double const p_b = 5000;
 
-    MPL::Property const& pressure_saturation =
-        MPL::SaturationVanGenuchten{"saturation", residual_liquid_saturation,
-                                    residual_gas_saturation, exponent, p_b};
+    MPL::Property const& pressure_saturation = MPL::SaturationVanGenuchten{
+        "saturation",      residual_liquid_saturation, residual_gas_saturation,
+        pressure_exponent, saturation_exponent,        p_b};
 
     MPL::VariableArray variable_array;
     ParameterLib::SpatialPosition const pos;

Tests/ProcessLib/RichardsMechanics/MicroporosityComputation.cpp

@@ -78,11 +78,12 @@ TEST(RichardsMechanics, computeMicroPorosity)
     // Saturation
     constexpr double residual_liquid_saturation = 0;
     constexpr double residual_gas_saturation = 0;
-    constexpr double exponent = 0.4;
+    constexpr double pressure_exponent = 0.4;
+    constexpr double saturation_exponent = 1. / (1. - pressure_exponent);
     constexpr double p_b = 300e6;
     MPL::Property const& saturation_micro = MPL::SaturationVanGenuchten{
         "saturation_micro", residual_liquid_saturation, residual_gas_saturation,
-        exponent, p_b};
+        pressure_exponent,  saturation_exponent,        p_b};
 
     // Swelling
     constexpr std::array swelling_pressures{1e7, 1e7, 1e7};