Merge branch 'improve_permeablility' into 'master'

[MPL/RelPermNonWettingPhaseVanGenuchtenMualem]  replace the Newton-Raphson method with the bisection method

See merge request ogs/ogs!3923

MaterialLib/MPL/Properties/RelativePermeability/CreateRelPermNonWettingPhaseVanGenuchtenMualem.cpp

@@ -44,6 +44,15 @@ std::unique_ptr<Property> createRelPermNonWettingPhaseVanGenuchtenMualem(
         //! \ogs_file_param{properties__property__RelativePermeabilityNonWettingPhaseVanGenuchtenMualem__min_relative_permeability}
         config.getConfigParameter<double>("min_relative_permeability");
 
+    if (min_relative_permeability <= 0.0 || min_relative_permeability > 1.0)
+    {
+        OGS_FATAL(
+            "The value for min_relative_permeability of "
+            "RelativePermeabilityNonWettingPhaseVanGenuchtenMualem is {:g}, "
+            "which falls outside of the range of (0, 1]",
+            min_relative_permeability);
+    }
+
     return std::make_unique<RelPermNonWettingPhaseVanGenuchtenMualem>(
         property_name, residual_liquid_saturation, residual_gas_saturation,
         exponent, min_relative_permeability);

MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.cpp

@@ -18,9 +18,17 @@
 #include "BaseLib/Error.h"
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Utils/CheckVanGenuchtenExponentRange.h"
+#include "MathLib/Nonlinear/Root1D.h"
 
 namespace MaterialPropertyLib
 {
+double computeVanGenuchtenMualemValue(const double S_L, const double S_L_r,
+                                      const double S_L_max, const double m)
+{
+    const double Se = (S_L - S_L_r) / (S_L_max - S_L_r);
+    return std::sqrt(1.0 - Se) * std::pow(1.0 - std::pow(Se, 1.0 / m), 2.0 * m);
+}
+
 RelPermNonWettingPhaseVanGenuchtenMualem::
     RelPermNonWettingPhaseVanGenuchtenMualem(std::string name,
                                              const double S_L_r,
@@ -47,7 +55,8 @@ PropertyDataType RelPermNonWettingPhaseVanGenuchtenMualem::value(
             variable_array[static_cast<int>(Variable::liquid_saturation)]),
         S_L_r_, S_L_max_);
 
-    const double krel = computeValue(S_L);
+    const double krel =
+        computeVanGenuchtenMualemValue(S_L, S_L_r_, S_L_max_, m_);
     return std::max(krel_min_, krel);
 }
 
@@ -70,20 +79,6 @@ PropertyDataType RelPermNonWettingPhaseVanGenuchtenMualem::dValue(
         return 0.0;
     }
 
-    return computeDerivative(S_L);
-}
-
-double RelPermNonWettingPhaseVanGenuchtenMualem::computeValue(
-    const double S_L) const
-{
-    const double Se = (S_L - S_L_r_) / (S_L_max_ - S_L_r_);
-    return std::sqrt(1.0 - Se) *
-           std::pow(1.0 - std::pow(Se, 1.0 / m_), 2.0 * m_);
-}
-
-double RelPermNonWettingPhaseVanGenuchtenMualem::computeDerivative(
-    const double S_L) const
-{
     if (std::fabs(S_L - S_L_max_) < std::numeric_limits<double>::epsilon())
     {
         return 0.0;
@@ -99,28 +94,23 @@ double RelPermNonWettingPhaseVanGenuchtenMualem::computeDerivative(
                 std::pow(val2, 2.0 * m_ - 1.0)) /
            (S_L_max_ - S_L_r_);
 }
+
 double RelPermNonWettingPhaseVanGenuchtenMualem::
     computeSaturationForMinimumRelativePermeability() const
 {
-    double S_for_k_rel_min = 0.5 * (S_L_r_ + S_L_max_);
-    const double tolerance = 1.e-16;
-    const double r0 = computeValue(S_for_k_rel_min) - krel_min_;
-    for (int iterations = 0; iterations < 1000; ++iterations)
+    auto f = [this](double S_L) -> double
     {
-        const double r = computeValue(S_for_k_rel_min) - krel_min_;
-        S_for_k_rel_min = std::clamp(S_for_k_rel_min, S_L_r_, S_L_max_);
-        S_for_k_rel_min -= r / computeDerivative(S_for_k_rel_min);
-
-        if (std::fabs(r / r0) < tolerance)
-        {
-            return S_for_k_rel_min;
-        }
-    }
+        return computeVanGenuchtenMualemValue(S_L, this->S_L_r_, this->S_L_max_,
+                                              this->m_) -
+               this->krel_min_;
+    };
+
+    auto root_finder =
+        MathLib::Nonlinear::makeRegulaFalsi<MathLib::Nonlinear::Pegasus>(
+            f, S_L_r_, S_L_max_);
+
+    root_finder.step(1000);
 
-    OGS_FATAL(
-        "The given minimum relative permeability, {:g}, is  not "
-        "associated with the saturation  in the range of '('{:f}, {:f}')'. "
-        "Please try another one in '(0, 1)', which should be close to zero",
-        krel_min_, S_L_r_, S_L_max_);
+    return root_finder.getResult();
 }
 }  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.h

@@ -33,6 +33,16 @@ class Medium;
  *       &S^L_{\mbox{max}}& \mbox{maximum saturation of wetting phase,}\\
  *       &m\, \in (0, 1) &    \mbox{ exponent.}\\
  *    \f}
+ *
+ *  The derivative of the relative permeability with respect to saturation is
+ *  computed as
+ *  \f[\frac{\mathrm{d} k_{rel}^n}{\mathrm{d}S^L}=
+ *  -(\dfrac{[1-S_e^{1/m}]^{2m}}{2\sqrt{1-S_e}}+2\sqrt{1-S_e}
+ *  {(1-S_e^{1/m})}^{2*m-1}
+ *   S_e^{1/m-1})/(S^L_\mbox{max}-S^L_r) \f]
+ *   As \f$S^L \to S^L_\mbox{max}\f$, or \f$S_e \to 1\f$,
+ *  \f$\dfrac{[1-S_e^{1/m}]^{2m}}{2\sqrt{1-S_e}}\f$ has a limit of zero.
+ *
  */
 class RelPermNonWettingPhaseVanGenuchtenMualem final : public Property
 {
@@ -74,6 +84,13 @@ public:
                             ParameterLib::SpatialPosition const& pos,
                             double const t, double const dt) const override;
 
+    /**
+     * Computes the saturation that gives the minimum relative permeability by
+     * using the Regula–Falsi Method.
+     * @return \f$ S^L\f$ that gives the minimum relative permeability.
+     */
+    double computeSaturationForMinimumRelativePermeability() const;
+
 private:
     const double S_L_r_;     ///< Residual saturation of wetting phase.
     const double S_L_max_;   ///< Maximum saturation of wetting phase.
@@ -81,31 +98,5 @@ private:
     const double krel_min_;  ///< Minimum relative permeability.
     const double
         S_L_for_krel_min_;  ///< Liquid saturation that gives \c krel_min_.
-
-    /**
-     *
-     * @param S_L Liquid saturation.
-     * @return \return \\return \f$k_{rel}^n \f$.
-     */
-    double computeValue(const double S_L) const;
-    /**
-     *  Computes
-     * \f[\frac{\mathrm{d} k_{rel}^n}{\mathrm{d}S^L}=
-     * -(\dfrac{[1-S_e^{1/m}]^{2m}}{2\sqrt{1-S_e}}+2\sqrt{1-S_e}
-     * {(1-S_e^{1/m})}^{2*m-1}
-     *  S_e^{1/m-1})/(S^L_\mbox{max}-S^L_r) \f]
-     *  As \f$S^L \to S^L_\mbox{max}\f$, or \f$S_e \to 1\f$,
-     * \f$\dfrac{[1-S_e^{1/m}]^{2m}}{2\sqrt{1-S_e}}\f$ has
-     * limit zero.
-     * @param S_L Liquid saturation.
-     * @return \return \f$ \frac{\mathrm{d} k_{rel}^n}{\mathrm{d}S^L} \f$.
-     */
-    double computeDerivative(const double S_L) const;
-    /**
-     * Computes the saturation that gives the minimum relative permeability by
-     * using the Newton-Raphson method.
-     * @return \f$ S^L\f$ that gives the minimum relative permeability.
-     */
-    double computeSaturationForMinimumRelativePermeability() const;
 };
 }  // namespace MaterialPropertyLib

Tests/MaterialLib/TestMPLRelPermNonWettingPhaseVanGenuchtenMualem.cpp

@@ -141,4 +141,12 @@ TEST(MaterialPropertyLib, RelPermNonWettingPhaseVanGenuchtenMualem)
         ASSERT_NEAR(dkrel_dS, factor * (k_rel_i_1 - k_rel_i_0) / perturbation,
                     6.0e-8);
     }
+
+    // Test the calculation of the liquid saturation for the minimum relative
+    // permeability.
+    MPL::RelPermNonWettingPhaseVanGenuchtenMualem k_non_wetting(
+        "dummy", 0.1, 0.4, 0.3288590604, 1.e-9);
+    ASSERT_NEAR(0.59999999552634464,
+                k_non_wetting.computeSaturationForMinimumRelativePermeability(),
+                1.0e-16);
 }