Merge pull request #2652 from Scinopode/mpl_brooks_corey_saturation

New MPL property: Brooks&Corey saturation curve

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

+Saturation model presented by Brooks&Corey (1964) \cite BrooksCorey1964.

Documentation/ProjectFile/properties/property/SaturationBrooksCorey/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/SaturationBrooksCorey/t_lambda.md

+The exponent of the Brooks&Corey \cite BrooksCorey1964 saturation function.

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

+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

+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"
\ No newline at end of file
+#include "CreateSaturationBrooksCorey.h"
\ No newline at end of file

MaterialLib/MPL/Properties/CreateSaturationBrooksCorey.cpp

+/**
+ * \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

+/**
+ * \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

+/**
+ * \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

+/**
+ * \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

Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp

+/**
+ * \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 <gtest/gtest.h>
+#include <sstream>
+
+#include "TestMPL.h"
+#include "Tests/TestTools.h"
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Properties/SaturationBrooksCorey.h"
+
+TEST(MaterialPropertyLib, SaturationBrooksCorey)
+{
+    const double ref_lambda = 2.5;
+    const double ref_residual_liquid_saturation = 0.12;
+    const double ref_residual_gas_saturation = 0.06;
+    const double ref_entry_pressure = 5678.54;
+
+    const double max_saturation = 1. - ref_residual_gas_saturation;
+
+    std::stringstream m;
+    m << "<medium>\n";
+    m << "<phases></phases>\n";
+    m << "<properties>\n";
+    m << "  <property>\n";
+    m << "    <name>saturation</name>\n";
+    m << "    <type>SaturationBrooksCorey</type>\n";
+    m << "    <residual_liquid_saturation>" << ref_residual_liquid_saturation
+      << "</residual_liquid_saturation>\n";
+    m << "    <residual_gas_saturation>" << ref_residual_gas_saturation
+      << "</residual_gas_saturation>\n";
+    m << "    <lambda>" << ref_lambda << "</lambda>\n";
+    m << "    <entry_pressure>" << ref_entry_pressure << "</entry_pressure>\n";
+    m << "  </property> \n";
+    m << "</properties>\n";
+    m << "</medium>\n";
+
+    auto const& medium = createTestMaterial(m.str());
+
+    MaterialPropertyLib::VariableArray variable_array;
+    ParameterLib::SpatialPosition const pos;
+    double const time = std::numeric_limits<double>::quiet_NaN();
+
+    variable_array[static_cast<int>(
+        MaterialPropertyLib::Variable::capillary_pressure)] = 0.0;
+
+    auto s_L = medium->property(MaterialPropertyLib::PropertyType::saturation)
+                   .template value<double>(variable_array, pos, time);
+    auto ds_L_dp_cap =
+        medium->property(MaterialPropertyLib::PropertyType::saturation)
+            .template dValue<double>(
+                variable_array,
+                MaterialPropertyLib::Variable::capillary_pressure,
+                pos,
+                time);
+
+    ASSERT_EQ(s_L, max_saturation);
+    ASSERT_EQ(ds_L_dp_cap, 0.);
+
+    for (double p_cap = 1.0; p_cap < 1.0e10; p_cap *= 1.5)
+    {
+        variable_array[static_cast<int>(
+            MaterialPropertyLib::Variable::capillary_pressure)] = p_cap;
+
+        s_L = medium->property(MaterialPropertyLib::PropertyType::saturation)
+                  .template value<double>(variable_array, pos, time);
+        ds_L_dp_cap =
+            medium->property(MaterialPropertyLib::PropertyType::saturation)
+                .template dValue<double>(
+                    variable_array,
+                    MaterialPropertyLib::Variable::capillary_pressure,
+                    pos,
+                    time);
+
+        const double s_eff =
+            std::pow(ref_entry_pressure / std::max(p_cap, ref_entry_pressure),
+                     ref_lambda);
+        const double s_ref =
+            s_eff * (max_saturation - ref_residual_liquid_saturation) +
+            ref_residual_liquid_saturation;
+        const double ds_eff_dpc =
+            (p_cap <= ref_entry_pressure) ? 0. : -ref_lambda / p_cap * s_ref;
+            const double ds_L_ds_eff = 1. / (max_saturation - ref_residual_liquid_saturation);
+            const double ds_L_dpc = ds_L_ds_eff * ds_eff_dpc;
+
+        ASSERT_NEAR(s_L, s_ref, 1.e-10);
+        ASSERT_NEAR(ds_L_dp_cap, ds_L_dpc, 1.e-10);
+    }
+}
\ No newline at end of file