diff --git a/MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.cpp b/MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.cpp
index bb824fb3d7eb3db9e9159149c74542af14b50ea4..3536e48b1a16cfd7a1a1ed0ced29b193923eeb47 100644
--- a/MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.cpp
+++ b/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
diff --git a/MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.h b/MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.h
index 8c35d7fa4b3ad4b305a055a75cd7fcca01f5fb8b..16b27d3f2e792dd9b63eb665f3f5664ae84b1214 100644
--- a/MaterialLib/MPL/Properties/RelativePermeability/RelPermNonWettingPhaseVanGenuchtenMualem.h
+++ b/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