diff --git a/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/c_SaturationBrooksCorey.md b/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/c_SaturationBrooksCorey.md
new file mode 100644
index 0000000000000000000000000000000000000000..e0be5d6c8b45ec729e063d38f48a355f68e68df9
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/c_SaturationBrooksCorey.md
@@ -0,0 +1 @@
+Saturation model presented by Brooks&Corey (1964) \cite BrooksCorey1964.
diff --git a/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_entry_pressure.md b/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_entry_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..e3d246498cefcd2402e2df304f0047534317a04a
--- /dev/null
+++ b/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).
diff --git a/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_lambda.md b/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_lambda.md
new file mode 100644
index 0000000000000000000000000000000000000000..6bb97b683fc7ba9dc1ee03065912d88f39128217
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_lambda.md
@@ -0,0 +1 @@
+The exponent of the Brooks&Corey \cite BrooksCorey1964 saturation function.
diff --git a/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_residual_gas_saturation.md b/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_residual_gas_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..d8631eb27f4627345443d3cb426330a84c325db1
--- /dev/null
+++ b/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.
diff --git a/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_residual_liquid_saturation.md b/Documentation/ProjectFile/properties/property/SaturationBrooksCorey/t_residual_liquid_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..95ce0ed4b319db3ecc0cfc87b42c611dcd572210
--- /dev/null
+++ b/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.
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 4d6eb09f375a007d9084b4b7985675d623c1a045..02646eddf901dc4a1e81856cb37c3287cdee0f23 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/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",
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 160408503fa0ddafaf007d94ec333a0dba18e730..646a6e2213ab21729599674913def917a1b292cb 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -16,4 +16,5 @@
#include "CreateExponentialProperty.h"
#include "CreateIdealGasLaw.h"
#include "CreateLinearProperty.h"
-#include "CreateParameterProperty.h"
\ No newline at end of file
+#include "CreateParameterProperty.h"
+#include "CreateSaturationBrooksCorey.h"
\ No newline at end of file
diff --git a/MaterialLib/MPL/Properties/CreateSaturationBrooksCorey.cpp b/MaterialLib/MPL/Properties/CreateSaturationBrooksCorey.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3bfcf096c69599d4b9924783a01e71ff294f8015
--- /dev/null
+++ b/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
diff --git a/MaterialLib/MPL/Properties/CreateSaturationBrooksCorey.h b/MaterialLib/MPL/Properties/CreateSaturationBrooksCorey.h
new file mode 100644
index 0000000000000000000000000000000000000000..42a8286b8d54914bfc7f9ba8012820cee7134156
--- /dev/null
+++ b/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
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index f8ca051af51b98fed3182f812161212dcb057abf..bd4b40595de10470c7f9ee174806a5f92839190d 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -17,3 +17,4 @@
#include "IdealGasLaw.h"
#include "LinearProperty.h"
#include "ParameterProperty.h"
+#include "SaturationBrooksCorey.h"
diff --git a/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp b/MaterialLib/MPL/Properties/SaturationBrooksCorey.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..54f87cacfc5165eda475ba99e602f24b52971a74
--- /dev/null
+++ b/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
diff --git a/MaterialLib/MPL/Properties/SaturationBrooksCorey.h b/MaterialLib/MPL/Properties/SaturationBrooksCorey.h
new file mode 100644
index 0000000000000000000000000000000000000000..2ee2c3c2c1d6cc1ec4efe543c905c16783cfff0a
--- /dev/null
+++ b/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
diff --git a/Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp b/Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..45700e5fdc8b9dccd15d6d9554a606a4121fa8b1
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLSaturationBrooksCorey.cpp
@@ -0,0 +1,97 @@
+/**
+ * \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