From b14255b7fa1baae35a7069eba7f17b9459b44d6c Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed, 27 Jan 2021 17:14:18 +0100
Subject: [PATCH] [MPL/RelPermUdell] Removed the non-wetting phase part

---
 .../CreateRelPermUdell.cpp                    | 21 +++----
 .../RelativePermeability/RelPermUdell.cpp     | 61 ++++++++-----------
 .../RelativePermeability/RelPermUdell.h       | 16 ++---
 3 files changed, 41 insertions(+), 57 deletions(-)

diff --git a/MaterialLib/MPL/Properties/RelativePermeability/CreateRelPermUdell.cpp b/MaterialLib/MPL/Properties/RelativePermeability/CreateRelPermUdell.cpp
index b2c8622c404..7adcc5a17cc 100644
--- a/MaterialLib/MPL/Properties/RelativePermeability/CreateRelPermUdell.cpp
+++ b/MaterialLib/MPL/Properties/RelativePermeability/CreateRelPermUdell.cpp
@@ -28,27 +28,22 @@ std::unique_ptr<RelPermUdell> createRelPermUdell(
     DBUG("Create RelPermUdell medium property {:s}.", property_name);
 
     auto const residual_liquid_saturation =
-        //! \ogs_file_param{properties__property__RelPermUdell__residual_liquid_saturation}
+        //! \ogs_file_param{properties__property__RelativePermeabilityUdell__residual_liquid_saturation}
         config.getConfigParameter<double>("residual_liquid_saturation");
     auto const residual_gas_saturation =
-        //! \ogs_file_param{properties__property__RelPermUdell__residual_gas_saturation}
+        //! \ogs_file_param{properties__property__RelativePermeabilityUdell__residual_gas_saturation}
         config.getConfigParameter<double>("residual_gas_saturation");
-    auto const min_relative_permeability_liquid =
-        //! \ogs_file_param{properties__property__RelPermUdell__min_relative_permeability_liquid}
-        config.getConfigParameter<double>("min_relative_permeability_liquid");
-    auto const min_relative_permeability_gas =
-        //! \ogs_file_param{properties__property__RelPermUdell__min_relative_permeability_gas}
-        config.getConfigParameter<double>("min_relative_permeability_gas");
+    auto const min_relative_permeability =
+        //! \ogs_file_param{properties__property__RelativePermeabilityUdell__min_relative_permeability}
+        config.getConfigParameter<double>("min_relative_permeability");
 
-    if ((min_relative_permeability_liquid < 0) ||
-        (min_relative_permeability_gas < 0))
+    if (min_relative_permeability < 0)
     {
-        OGS_FATAL("Minimal relative permeabilities must be non-negative.");
+        OGS_FATAL("Minimal relative permeability must be non-negative.");
     }
 
     return std::make_unique<RelPermUdell>(
         std::move(property_name), residual_liquid_saturation,
-        residual_gas_saturation, min_relative_permeability_liquid,
-        min_relative_permeability_gas);
+        residual_gas_saturation, min_relative_permeability);
 }
 }  // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.cpp b/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.cpp
index 92a2ead88b9..3f046a6eef8 100644
--- a/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.cpp
+++ b/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.cpp
@@ -23,12 +23,10 @@ namespace MaterialPropertyLib
 RelPermUdell::RelPermUdell(std::string name,
                            const double residual_liquid_saturation,
                            const double residual_gas_saturation,
-                           const double min_relative_permeability_liquid,
-                           const double min_relative_permeability_gas)
+                           const double min_relative_permeability)
     : residual_liquid_saturation_(residual_liquid_saturation),
       residual_gas_saturation_(residual_gas_saturation),
-      min_relative_permeability_liquid_(min_relative_permeability_liquid),
-      min_relative_permeability_gas_(min_relative_permeability_gas)
+      min_relative_permeability_(min_relative_permeability)
 {
     name_ = std::move(name);
 }
@@ -38,69 +36,58 @@ PropertyDataType RelPermUdell::value(
     ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
     double const /*dt*/) const
 {
-    const double s_L = std::get<double>(
+    const double S_L = std::get<double>(
         variable_array[static_cast<int>(Variable::liquid_saturation)]);
 
-    if (std::isnan(s_L))
+    if (std::isnan(S_L))
     {
         OGS_FATAL("Liquid saturation not set in RelPermUdell::value().");
     }
 
-    auto const s_L_res = residual_liquid_saturation_;
-    auto const s_L_max = 1. - residual_gas_saturation_;
-    auto const k_rel_min_LR = min_relative_permeability_liquid_;
-    auto const k_rel_min_GR = min_relative_permeability_gas_;
+    auto const S_L_res = residual_liquid_saturation_;
+    auto const S_L_max = 1. - residual_gas_saturation_;
 
-    auto const s = (s_L - s_L_res) / (s_L_max - s_L_res);
+    auto const S_e = (S_L - S_L_res) / (S_L_max - S_L_res);
 
-    if (s >= 1.0)
+    if (S_e >= 1.0)
     {
         // fully saturated medium
-        return Eigen::Vector2d{1.0, k_rel_min_GR};
+        return 1.0;
     }
-    if (s <= 0.0)
+    if (S_e <= 0.0)
     {
         // dry medium
-        return Eigen::Vector2d{k_rel_min_LR, 1.0};
+        return min_relative_permeability_;
     }
 
-    auto const k_rel_LR = s * s * s;
-    auto const k_rel_GR = (1. - s) * (1. - s) * (1. - s);
-
-    return Eigen::Vector2d{std::max(k_rel_LR, k_rel_min_LR),
-                           std::max(k_rel_GR, k_rel_min_GR)};
+    return std::max(min_relative_permeability_, S_e * S_e * S_e);
 }
 PropertyDataType RelPermUdell::dValue(
-    VariableArray const& variable_array, Variable const primary_variable,
+    VariableArray const& variable_array, Variable const variable,
     ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
     double const /*dt*/) const
 {
-    (void)primary_variable;
-    assert((primary_variable == Variable::liquid_saturation) &&
+    (void)variable;
+    assert((variable == Variable::liquid_saturation) &&
            "RelPermUdell::dValue is implemented for "
            " derivatives with respect to liquid saturation only.");
 
-    const double s_L = std::get<double>(
+    const double S_L = std::get<double>(
         variable_array[static_cast<int>(Variable::liquid_saturation)]);
 
-    auto const s_L_res = residual_liquid_saturation_;
-    auto const s_L_max = 1. - residual_gas_saturation_;
-    auto const s = (s_L - s_L_res) / (s_L_max - s_L_res);
+    auto const S_L_res = residual_liquid_saturation_;
+    auto const S_L_max = 1. - residual_gas_saturation_;
+    auto const S_e = (S_L - S_L_res) / (S_L_max - S_L_res);
 
-    if ((s < 0.) || (s > 1.))
+    if ((S_e < 0.) || (S_e > 1.))
     {
-        return Eigen::Vector2d{0., 0.};
+        return 0.;
     }
 
-    auto const d_se_d_sL = 1. / (s_L_max - s_L_res);
-
-    auto const dk_rel_LRdse = 3. * s * s;
-    auto const dk_rel_LRdsL = dk_rel_LRdse * d_se_d_sL;
-
-    auto const dk_rel_GRdse = -3. * (1. - s) * (1. - s);
-    auto const dk_rel_GRdsL = dk_rel_GRdse * d_se_d_sL;
+    auto const dS_e_dS_L = 1. / (S_L_max - S_L_res);
 
-    return Eigen::Vector2d{dk_rel_LRdsL, dk_rel_GRdsL};
+    auto const dk_rel_LR_dS_e = 3. * S_e * S_e;
+    return dk_rel_LR_dS_e * dS_e_dS_L;
 }
 
 }  // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.h b/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.h
index 845c0a857e5..3eb7b2ec13f 100644
--- a/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.h
+++ b/MaterialLib/MPL/Properties/RelativePermeability/RelPermUdell.h
@@ -31,6 +31,10 @@ class Component;
  *          - \f$S^{\text{eff}}_{\alpha}\f$ is the effective saturation of
  * phase \f$\alpha\f$
  *
+ * This class handles the wetting (liquid) phase portion of this relative
+ * permeability property, i,e with \f$\alpha=L\f$ for the  relative permeability
+ * function.
+ *
  * \details This property must be a medium property, it
  * computes the permeability reduction due to saturation as function of
  * capillary pressure.
@@ -40,22 +44,20 @@ class RelPermUdell final : public Property
 private:
     const double residual_liquid_saturation_;
     const double residual_gas_saturation_;
-    const double min_relative_permeability_liquid_;
-    const double min_relative_permeability_gas_;
+    const double min_relative_permeability_;
 
 public:
     RelPermUdell(std::string name, const double residual_liquid_saturation,
                  const double residual_gas_saturation,
-                 const double min_relative_permeability_liquid,
-                 const double min_relative_permeability_gas);
+                 const double min_relative_permeability);
 
     void checkScale() const override
     {
         if (!std::holds_alternative<Medium*>(scale_))
         {
             OGS_FATAL(
-                "The property 'RelPermUdell' is implemented on the "
-                "'media' scale only.");
+                "The property 'RelativePermeabilityUdell' is implemented on "
+                "the 'media' scale only.");
         }
     }
 
@@ -63,7 +65,7 @@ public:
                            ParameterLib::SpatialPosition const& pos,
                            double const t, double const dt) const override;
     PropertyDataType dValue(VariableArray const& variable_array,
-                            Variable const primary_variable,
+                            Variable const variable,
                             ParameterLib::SpatialPosition const& pos,
                             double const t, double const dt) const override;
 };
-- 
GitLab