Brooks-Corey saturation model added to MPL properties

MaterialLib/MPL/CreateProperty.cpp

@@ -58,6 +58,12 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
     {
         return createIdealGasLaw(config);
     }
+
+    if (boost::iequals(property_type, "SaturationBrooksCorey"))
+    {
+        return createSaturationBrooksCorey(config);
+    }
+
     /* TODO Additional properties go here, for example:
     if (boost::iequals(property_type, "BilinearTemperaturePressure"))
     {

MaterialLib/MPL/Properties/Properties.h

@@ -17,3 +17,4 @@
 #include "IdealGasLaw.h"
 #include "LinearProperty.h"
 #include "ParameterProperty.h"
+#include "SaturationBrooksCorey.h"

MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp

+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date   27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, 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 "SaturationBrooksCorey.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MathLib/MathTools.h"
+
+namespace MaterialPropertyLib
+{
+SaturationBrooksCorey::SaturationBrooksCorey(
+    const double residual_liquid_saturation,
+    const double residual_gas_saturation,
+    const double exponent,
+    const double entry_pressure)
+    : _residual_liquid_saturation(residual_liquid_saturation),
+      _residual_gas_saturation(residual_gas_saturation),
+      _exponent(exponent),
+      _entry_pressure(entry_pressure){};
+
+PropertyDataType SaturationBrooksCorey::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/,
+    double const /*t*/) const
+{
+    const double p_cap = std::max(
+        std::numeric_limits<double>::min(),
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::capillary_pressure)]));
+
+    const double s_L_res = _residual_liquid_saturation;
+    const double s_L_max = 1.0 - _residual_gas_saturation;
+    const double lambda = _exponent;
+    const double p_b = _entry_pressure;
+
+    if (p_cap <= p_b)
+        return s_L_max;
+
+    const double s_eff = std::pow(p_b / p_cap, lambda);
+    const double s = s_eff * (s_L_max - s_L_res) + s_L_res;
+
+    return MathLib::limitValueInInterval(s, s_L_res + _minor_offset,
+                                         s_L_max - _minor_offset);
+}
+
+PropertyDataType SaturationBrooksCorey::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) &&
+           "SaturationBrooksCorey::dValue is implemented for "
+           " derivatives with respect to capillary pressure only.");
+
+    const double p_b = _entry_pressure;
+    const double p_cap = std::max(
+        p_b,
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::capillary_pressure)]));
+
+    auto const s_L = _medium->property(PropertyType::saturation)
+                         .template value<double>(variable_array, pos, t);
+
+    const double lambda = _exponent;
+
+    return -lambda / p_cap * s_L;
+}
+
+PropertyDataType SaturationBrooksCorey::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) &&
+           "SaturationBrooksCorey::d2Value is implemented for "
+           " derivatives with respect to capillary pressure only.");
+
+    const double p_b = _entry_pressure;
+    const double p_cap = std::max(
+        p_b,
+        std::get<double>(
+            variable_array[static_cast<int>(Variable::capillary_pressure)]));
+
+    const double s_L_res = _residual_liquid_saturation;
+    const double s_L_max = 1.0 - _residual_gas_saturation;
+
+    const double lambda = _exponent;
+
+    return lambda * (lambda + 1) * std::pow(p_b / p_cap, lambda) /
+           (p_cap * p_cap) * (s_L_max - s_L_res);
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/SaturationBrooksCorey.h

+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date   27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2019, 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 <limits>
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+/**
+ * \brief A well known soil characteristics function
+ * \details This property must be a medium property, it
+ * computes the saturation of the wetting phase as function
+ * of capillary pressure.
+ */
+class SaturationBrooksCorey final : public Property
+{
+private:
+    Medium* _medium = nullptr;
+    const double _residual_liquid_saturation;
+    const double _residual_gas_saturation;
+    const double _exponent;
+    const double _entry_pressure;
+
+public:
+    SaturationBrooksCorey(const double residual_liquid_saturation,
+                          const double residual_gas_saturation,
+                          const double exponent,
+                          const double entry_pressure);
+
+    void setScale(
+        std::variant<Medium*, Phase*, Component*> scale_pointer) override
+    {
+        if (!std::holds_alternative<Medium*>(scale_pointer))
+        {
+            OGS_FATAL(
+                "The property 'SaturationBrooksCorey' is implemented on the "
+                "'media' scale only.");
+        }
+        _medium = std::get<Medium*>(scale_pointer);
+    }
+
+    /// 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;
+    PropertyDataType d2Value(VariableArray const& variable_array,
+                             Variable const variable1,
+                             Variable const variable2,
+                             ParameterLib::SpatialPosition const& /*pos*/,
+                             double const /*t*/) const override;
+};
+
+inline std::unique_ptr<SaturationBrooksCorey> createSaturationBrooksCorey(
+    BaseLib::ConfigTree const& config)
+{
+    // check is reading the parameter, not peeking it...
+    //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey}
+    // config.checkConfigParameter("type", "SaturationBrooksCorey");
+    DBUG("Create SaturationBrooksCorey medium property");
+
+    auto const residual_liquid_saturation =
+        //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__residual_liquid_saturation}
+        config.getConfigParameter<double>("residual_liquid_saturation");
+    auto const residual_gas_saturation =
+        //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__residual_gas_saturation}
+        config.getConfigParameter<double>("residual_gas_saturation");
+    auto const exponent =
+        //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__lambda}
+        config.getConfigParameter<double>("lambda");
+    auto const entry_pressure =
+        //! \ogs_file_param{prj__media__medium__properties__property__SaturationBrooksCorey__entry_pressure}
+        config.getConfigParameter<double>("entry_pressure");
+
+    return std::make_unique<SaturationBrooksCorey>(residual_liquid_saturation,
+                                                   residual_gas_saturation,
+                                                   exponent, entry_pressure);
+}
+
+}  // namespace MaterialPropertyLib