New MPL property: Brooks&Corey saturation curve

Documentation/ProjectFile/properties/property/SaturationBrooksCorey/c_SaturationBrooksCorey.md

@@ -0,0 +1 @@
+Saturation model presented by Brooks&Corey (1964) \cite BrooksCorey1964.

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

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

Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_lambda.md

@@ -0,0 +1 @@
+The exponent of the Brooks&Corey \cite BrooksCorey1964 saturation function.

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

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

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

@@ -0,0 +1 @@
+The smallest degree of saturation of the liquid (or wetting) phase in the medium.

MaterialLib/MPL/CreateProperty.cpp

@@ -18,6 +18,7 @@
 #include "BaseLib/ConfigTree.h"
 
 #include "Properties/CreateProperties.h"
+
 #include "Properties/Properties.h"
 
 #include "Component.h"
@@ -58,12 +59,11 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
     {
         return createIdealGasLaw(config);
     }
-    /* TODO Additional properties go here, for example:
-    if (boost::iequals(property_type, "BilinearTemperaturePressure"))
+
+    if (boost::iequals(property_type, "SaturationBrooksCorey"))
     {
-        return createBilinearTemperaturePressure(config, material_type);
+        return createSaturationBrooksCorey(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",

MaterialLib/MPL/Properties/CreateProperties.h

@@ -16,4 +16,5 @@
 #include "CreateExponentialProperty.h"
 #include "CreateIdealGasLaw.h"
 #include "CreateLinearProperty.h"
-#include "CreateParameterProperty.h"
+#include "CreateParameterProperty.h"
+#include "CreateSaturationBrooksCorey.h"

MaterialLib/MPL/Properties/CreateSaturationBrooksCorey.cpp

@@ -0,0 +1,42 @@
+/**
+ * \file
+ * \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 "BaseLib/ConfigTree.h"
+#include "SaturationBrooksCorey.h"
+
+namespace MaterialPropertyLib
+{
+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

MaterialLib/MPL/Properties/CreateSaturationBrooksCorey.h

@@ -0,0 +1,27 @@
+/**
+ * \file
+ * \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
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class SaturationBrooksCorey;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationBrooksCorey> createSaturationBrooksCorey(
+    BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

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

@@ -0,0 +1,109 @@
+/**
+ * \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::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);
+    return s_eff * (s_L_max - s_L_res) + s_L_res;
+}
+
+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 s_L_res = _residual_liquid_saturation;
+    const double s_L_max = 1.0 - _residual_gas_saturation;
+    const double p_b = _entry_pressure;
+    const double p_cap = std::get<double>(
+        variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+    if (p_cap <= p_b)
+        return 0.;
+
+    auto const s_L = _medium->property(PropertyType::saturation)
+                         .template value<double>(variable_array, pos, t);
+
+    const double lambda = _exponent;
+    const double ds_L_d_s_eff = 1. / (s_L_max - s_L_res);
+
+    return -lambda / p_cap * s_L * ds_L_d_s_eff;
+}
+
+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

@@ -0,0 +1,72 @@
+/**
+ * \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 "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;
+};
+}  // namespace MaterialPropertyLib