Merge pull request #2788 from endJunction/MPLVanGenuchtenSaturationPermeability

MPL van Genuchten saturation and relative permeability models.

Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/c_RelativePermeabilityVanGenuchten.md

+\copydoc MaterialPropertyLib::RelativePermeabilityVanGenuchten

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

+The exponent of the van Genuchten relative permeability function.

Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_minimum_relative_permeability_liquid.md

+The minimum relative permeability of the liquid (or wetting) phase.

Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_gas_saturation.md

+The smallest degree of saturation of the gas (or non-wetting) phase in the medium.

Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_liquid_saturation.md

+The smallest degree of saturation of the liquid (or wetting) phase in the medium.

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

+\copydoc MaterialPropertyLib::SaturationVanGenuchten

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

+The required pressure for a non-wetting fluid to enter a porous medium (the capillary pressure at full saturation).

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

+The exponent of the van Genuchten saturation function.

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

+The smallest degree of saturation of the gas (or non-wetting) phase in the medium.

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

+The smallest degree of saturation of the liquid (or wetting) phase in the medium.

MaterialLib/MPL/CreateProperty.cpp

@@ -79,6 +79,16 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
         return createRelPermLiakopoulos(config);
     }
 
+    if (boost::iequals(property_type, "SaturationVanGenuchten"))
+    {
+        return createSaturationVanGenuchten(config);
+    }
+
+    if (boost::iequals(property_type, "RelativePermeabilityVanGenuchten"))
+    {
+        return createRelPermVanGenuchten(config);
+    }
+
     // If none of the above property types are found, OGS throws an error.
     OGS_FATAL("The specified component property type '%s' was not recognized",
               property_type.c_str());

MaterialLib/MPL/Properties/CreateProperties.h

@@ -18,5 +18,7 @@
 #include "CreateParameterProperty.h"
 #include "CreateRelPermBrooksCorey.h"
 #include "CreateRelPermLiakopoulos.h"
+#include "CreateRelPermVanGenuchten.h"
 #include "CreateSaturationBrooksCorey.h"
-#include "CreateSaturationLiakopoulos.h"
+#include "CreateSaturationLiakopoulos.h"
\ No newline at end of file
+#include "CreateSaturationVanGenuchten.h"

MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "RelPermVanGenuchten.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermVanGenuchten> createRelPermVanGenuchten(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "RelativePermeabilityVanGenuchten");
+    DBUG("Create RelativePermeabilityVanGenuchten medium property");
+
+    auto const residual_liquid_saturation =
+        //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__residual_liquid_saturation}
+        config.getConfigParameter<double>("residual_liquid_saturation");
+    auto const residual_gas_saturation =
+        //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__residual_gas_saturation}
+        config.getConfigParameter<double>("residual_gas_saturation");
+    auto const min_relative_permeability_liquid =
+        //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__minimum_relative_permeability_liquid}
+        config.getConfigParameter<double>(
+            "minimum_relative_permeability_liquid");
+    auto const exponent =
+        //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__exponent}
+        config.getConfigParameter<double>("exponent");
+    if (exponent <= 0. || exponent >= 1.)
+    {
+        OGS_FATAL("Exponent must be in the (0, 1) range.");
+    }
+
+    return std::make_unique<RelPermVanGenuchten>(
+        residual_liquid_saturation,
+        residual_gas_saturation,
+        min_relative_permeability_liquid,
+        exponent);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class RelPermVanGenuchten;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermVanGenuchten> createRelPermVanGenuchten(
+    BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "SaturationVanGenuchten.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "SaturationVanGenuchten");
+
+    DBUG("Create SaturationVanGenuchten medium property");
+
+    auto const residual_liquid_saturation =
+        //! \ogs_file_param{properties__property__SaturationVanGenuchten__residual_liquid_saturation}
+        config.getConfigParameter<double>("residual_liquid_saturation");
+    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");
+    auto const entry_pressure =
+        //! \ogs_file_param{properties__property__SaturationVanGenuchten__entry_pressure}
+        config.getConfigParameter<double>("entry_pressure");
+
+    return std::make_unique<SaturationVanGenuchten>(residual_liquid_saturation,
+                                                    residual_gas_saturation,
+                                                    exponent, entry_pressure);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class SaturationVanGenuchten;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten(
+    BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/Properties.h

@@ -18,5 +18,7 @@
 #include "ParameterProperty.h"
 #include "RelPermBrooksCorey.h"
 #include "RelPermLiakopoulos.h"
+#include "RelPermVanGenuchten.h"
 #include "SaturationBrooksCorey.h"
 #include "SaturationLiakopoulos.h"
+#include "SaturationVanGenuchten.h"

MaterialLib/MPL/Properties/RelPermVanGenuchten.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "RelPermVanGenuchten.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+RelPermVanGenuchten::RelPermVanGenuchten(
+    double const residual_liquid_saturation,
+    double const residual_gas_saturation,
+    double const min_relative_permeability_liquid,
+    double const exponent)
+    : _S_L_res(residual_liquid_saturation),
+      _S_L_max(1. - residual_gas_saturation),
+      _k_rel_min(min_relative_permeability_liquid),
+      _m(exponent)
+{
+    if (!(_m > 0 && _m < 1))
+    {
+        OGS_FATAL(
+            "The exponent value m = %g of van Genuchten relative permeability "
+            "model, is out of its range of (0, 1)",
+            _m);
+    }
+}
+
+PropertyDataType RelPermVanGenuchten::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/,
+    double const /*t*/) const
+{
+    double const S_L = std::clamp(
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::liquid_saturation)]),
+        _S_L_res, _S_L_max);
+
+    double const S_eff = (S_L - _S_L_res) / (_S_L_max - _S_L_res);
+    double const v = 1. - std::pow(1. - std::pow(S_eff, 1. / _m), _m);
+    double const k_rel = std::sqrt(S_eff) * v * v;
+    return std::max(_k_rel_min, k_rel);
+}
+
+PropertyDataType RelPermVanGenuchten::dValue(
+    VariableArray const& variable_array, Variable const primary_variable,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/) const
+{
+    (void)primary_variable;
+    assert((primary_variable == Variable::liquid_saturation) &&
+           "RelativePermeabilityVanGenuchten::dValue is implemented for "
+           "derivatives with respect to liquid saturation only.");
+    double const S_L = std::clamp(
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::liquid_saturation)]),
+        _S_L_res, _S_L_max);
+
+    double const S_eff = (S_L - _S_L_res) / (_S_L_max - _S_L_res);
+    if (S_eff <= 0)  // prevent division by zero
+    {
+        return 0;
+    }
+
+    if (S_eff >= 1)  // prevent taking root of zero
+    {
+        return 0;
+    }
+
+    double const S_eff_to_1_over_m = std::pow(S_eff, 1. / _m);
+    double const v = 1. - std::pow(1. - S_eff_to_1_over_m, _m);
+    double const sqrt_S_eff = std::sqrt(S_eff);
+    double const k_rel = sqrt_S_eff * v * v;
+
+    if (k_rel < _k_rel_min)
+    {
+        return 0;
+    }
+
+    return (0.5 * v * v / sqrt_S_eff +
+            2. * sqrt_S_eff * v * std::pow(1. - S_eff_to_1_over_m, _m - 1.) *
+                S_eff_to_1_over_m / S_eff) /
+           (_S_L_max - _S_L_res);
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/RelPermVanGenuchten.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+#pragma once
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+
+/// Van Genuchten relative permeability function.
+///
+/// This property must be a medium property, it computes the permeability
+/// reduction due to saturation as function of capillary pressure.
+class RelPermVanGenuchten final : public Property
+{
+private:
+    double const _S_L_res;
+    double const _S_L_max;
+    double const _k_rel_min;
+    double const _m;
+    Medium* _medium = nullptr;
+
+public:
+    RelPermVanGenuchten(double const residual_liquid_saturation,
+                        double const residual_gas_saturation,
+                        double const min_relative_permeability_liquid,
+                        double const exponent);
+    /// This method assigns a pointer to the material object that is the owner
+    /// of this property
+    void setScale(
+        std::variant<Medium*, Phase*, Component*> scale_pointer) override
+    {
+        if (std::holds_alternative<Medium*>(scale_pointer))
+        {
+            _medium = std::get<Medium*>(scale_pointer);
+        }
+        else
+        {
+            OGS_FATAL(
+                "The property 'RelativePermeabilityVanGenuchten' is "
+                "implemented on the 'media' scale only.");
+        }
+    }
+
+    /// Those methods override the base class implementations and
+    /// actually compute and set the property _values and _dValues.
+    PropertyDataType value(VariableArray const& variable_array,
+                           ParameterLib::SpatialPosition const& pos,
+                           double const t) const override;
+    PropertyDataType dValue(VariableArray const& variable_array,
+                            Variable const variable,
+                            ParameterLib::SpatialPosition const& pos,
+                            double const t) const override;
+};
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/SaturationVanGenuchten.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "SaturationVanGenuchten.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+SaturationVanGenuchten::SaturationVanGenuchten(
+    double const residual_liquid_saturation,
+    double const residual_gas_saturation,
+    double const exponent,
+    double const entry_pressure)
+    : _S_L_res(residual_liquid_saturation),
+      _S_L_max(1. - residual_gas_saturation),
+      _m(exponent),
+      _p_b(entry_pressure)
+{
+    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)",
+            _m);
+    }
+}
+
+PropertyDataType SaturationVanGenuchten::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/,
+    double const /*t*/) const
+{
+    const double p_cap = std::get<double>(
+        variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+    if (p_cap <= 0)
+    {
+        return _S_L_max;
+    }
+
+    double const p = p_cap / _p_b;
+    double const n = 1. / (1. - _m);
+    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;
+    return std::clamp(S, _S_L_res, _S_L_max);
+}
+
+PropertyDataType SaturationVanGenuchten::dValue(
+    VariableArray const& variable_array, Variable const primary_variable,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/) const
+{
+    (void)primary_variable;
+    assert((primary_variable == Variable::capillary_pressure) &&
+           "SaturationVanGenuchten::dvalue is implemented for derivatives with "
+           "respect to capillary pressure only.");
+
+    const double p_cap = std::get<double>(
+        variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+    if (p_cap <= 0)
+    {
+        return 0;
+    }
+
+    double const p = p_cap / _p_b;
+    double const n = 1. / (1. - _m);
+    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;
+
+    if (S < _S_L_res || S > _S_L_max)
+    {
+        return 0;
+    }
+
+    double const dS_eff_dp_cap = -_m * std::pow(p_cap / _p_b, n - 1) *
+                                 std::pow(1 + p_to_n, -1. - _m) /
+                                 (_p_b * (1. - _m));
+    return dS_eff_dp_cap * (_S_L_max - _S_L_res);
+}
+
+PropertyDataType SaturationVanGenuchten::d2Value(
+    VariableArray const& variable_array, Variable const primary_variable1,
+    Variable const primary_variable2,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/) const
+{
+    (void)primary_variable1;
+    (void)primary_variable2;
+    assert((primary_variable1 == Variable::capillary_pressure) &&
+           (primary_variable2 == Variable::capillary_pressure) &&
+           "SaturationVanGenuchten::d2Value is implemented for "
+           " derivatives with respect to capillary pressure only.");
+
+    const double p_cap = std::get<double>(
+        variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+    if (p_cap <= 0)
+    {
+        return 0;
+    }
+
+    double const p = p_cap / _p_b;
+    double const n = 1. / (1. - _m);
+    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;
+
+    if (S < _S_L_res || S > _S_L_max)
+    {
+        return 0;
+    }
+
+    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.));
+    return d2S_eff_dp_cap2 * (_S_L_max - _S_L_res);
+}
+}  // namespace MaterialPropertyLib