diff --git a/MaterialLib/CMakeLists.txt b/MaterialLib/CMakeLists.txt
index 60c46c969befbd935b7363f3493e47e78e6107bb..123a4514098c4718a943572ed5176e64aa4dad7c 100644
--- a/MaterialLib/CMakeLists.txt
+++ b/MaterialLib/CMakeLists.txt
@@ -11,6 +11,7 @@ append_source_files(SOURCES Fluid/Viscosity)
append_source_files(SOURCES PorousMedium/Porosity)
append_source_files(SOURCES PorousMedium/Storage)
append_source_files(SOURCES PorousMedium/Permeability)
+append_source_files(SOURCES PorousMedium/UnsaturatedProperty/CapillaryPressure)
add_library(MaterialLib ${SOURCES} )
target_link_libraries(MaterialLib
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9d8502961542f3178945bfe619ecff755d65e23a
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.cpp
@@ -0,0 +1,54 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file BrookCoreyCapillaryPressure.cpp
+ *
+ * Created on November 1, 2016, 9:45 AM
+ */
+
+#include "BrookCoreyCapillaryPressureSaturation.h"
+
+#include <cmath>
+
+#include "MathLib/MathTools.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+double BrookCoreyCapillaryPressureSaturation::getCapillaryPressure(
+ const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double Se = (S - _Sr) / (_Smax - _Sr);
+ const double pc = _pb * std::pow(Se, -1.0 / _mm);
+ return MathLib::limitValueInInterval(pc, _minor_offset, _Pc_max);
+}
+
+double BrookCoreyCapillaryPressureSaturation::getSaturation(
+ const double capillary_pressure) const
+{
+ const double pc =
+ (capillary_pressure < 0.0) ? _minor_offset : capillary_pressure;
+ const double Se = std::pow(pc / _pb, -_mm);
+ const double S = Se * (_Smax - _Sr) + _Sr;
+ return MathLib::limitValueInInterval(S, _Sr + _minor_offset,
+ _Smax - _minor_offset);
+}
+
+double BrookCoreyCapillaryPressureSaturation::getdPcdS(
+ const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double val = std::pow(((S - _Sr) / (_Smax - _Sr)), -1.0 / _mm);
+ return (_pb * val) / (_mm * (_Sr - S));
+}
+
+} // end namespace
+} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h
new file mode 100644
index 0000000000000000000000000000000000000000..aeb24f13ab5d866edf8258e164794a3cc6b6b183
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h
@@ -0,0 +1,92 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file BrookCoreyCapillaryPressureSaturation.h
+ *
+ * Created on November 1, 2016, 9:45 AM
+ */
+
+#ifndef OGS_BROOK_COREY_CAPILLARY_PRESSURE_SATURATION_H
+#define OGS_BROOK_COREY_CAPILLARY_PRESSURE_SATURATION_H
+
+#include <limits> // std::numeric_limits
+
+#include "CapillaryPressureSaturation.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+/**
+ * \brief Brook-Corey capillary pressure saturation model
+ *
+ * \f[p_c=p_b S_e^{-1/m}\f]
+ * with
+ * \f[S_e=\dfrac{S-S_r}{S_{\mbox{max}}-S_r}\f]
+ * where
+ * \f{eqnarray*}{
+ * &p_b& \mbox{ entry pressure,}\\
+ * &S_r& \mbox{ residual saturation,}\\
+ * &S_{\mbox{max}}& \mbox{ maximum saturation,}\\
+ * &m(>=1) & \mbox{ exponent.}\\
+ * \f}
+ *
+ * Note:
+ * \f[m=1/(1-n)\f]
+ */
+class BrookCoreyCapillaryPressureSaturation final
+ : public CapillaryPressureSaturation
+{
+public:
+ /**
+ * @param pb Entry pressure, \f$ p_b \f$
+ * @param Sr Residual saturation, \f$ S_r \f$
+ * @param Smax Maximum saturation, \f$ S_{\mbox{max}} \f$
+ * @param m Exponent, \f$ m \f$
+ * @param Pc_max Maximum capillary pressure, \f$ P_c^{\mbox{max}}\f$
+ */
+ BrookCoreyCapillaryPressureSaturation(const double pb, const double Sr,
+ const double Smax, const double m,
+ const double Pc_max)
+ : _pb(pb), _Sr(Sr), _Smax(Smax), _mm(m), _Pc_max(Pc_max)
+ {
+ }
+
+ /// Get model name.
+ std::string getName() const override
+ {
+ return "Brook-Corey capillary pressure saturation model.";
+ }
+
+ /// Get capillary pressure.
+ double getCapillaryPressure(const double saturation) const override;
+
+ /// Get capillary pressure.
+ double getSaturation(const double capillary_pressure) const override;
+
+ /// Get the derivative of the capillary pressure with respect to saturation
+ double getdPcdS(const double saturation) const override;
+
+private:
+ const double _pb; ///< Entry pressure.
+ const double _Sr; ///< Residual saturation.
+ const double _Smax; ///< Maximum saturation.
+ const double _mm; ///< Exponent (<=1.0), n=1/(1-mm).
+ const double _Pc_max; ///< Maximum capillaray pressure
+
+ /** A small number for an offset:
+ * 1. to set the bound of S, the saturation, such that
+ * S in [_Sr+_minor_offset, _Smax-_minor_offset]
+ * 2. to set the bound of Pc, the capillary pressure, such that
+ * Pc in [_minor_offset, _Pc_max]
+ */
+ const double _minor_offset = std::numeric_limits<double>::epsilon();
+};
+
+} // end namespace
+} // end namespace
+#endif /* OGS_BROOK_COREY_CAPILLARY_PRESSURE_SATURATION_H */
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h
new file mode 100644
index 0000000000000000000000000000000000000000..3319c41a5959c24c4e40f33cd122da9961c1ef16
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h
@@ -0,0 +1,43 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file: CapillaryPressureSaturation.h
+ *
+ */
+
+#ifndef OGS_CAPILLARY_PRESSURE_SATURATION_H
+#define OGS_CAPILLARY_PRESSURE_SATURATION_H
+
+#include <string>
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+/// Base class of capillary pressure models
+class CapillaryPressureSaturation
+{
+public:
+ virtual ~CapillaryPressureSaturation() = default;
+
+ /// Get model name.
+ virtual std::string getName() const = 0;
+
+ /// Get capillary pressure.
+ virtual double getCapillaryPressure(const double saturation) const = 0;
+
+ /// Get capillary pressure.
+ virtual double getSaturation(const double capillary_ressure) const = 0;
+
+ /// Get the derivative of the capillary pressure with respect to saturation
+ virtual double getdPcdS(const double saturation) const = 0;
+};
+
+} // end namespace
+} // end namespace
+
+#endif /* OGS_CAPILLARY_PRESSURE_SATURATION_H */
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d6e317e678746e6fd9b07dbc3f9422ce82619928
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp
@@ -0,0 +1,120 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateCapillaryPressureModel.cpp
+ *
+ * Created on November 1, 2016, 10:06 AM
+ */
+
+#include "CreateCapillaryPressureModel.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "BaseLib/Error.h"
+
+#include "CapillaryPressureSaturation.h"
+#include "BrookCoreyCapillaryPressureSaturation.h"
+#include "VanGenuchtenCapillaryPressureSaturation.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+/**
+ \param config ConfigTree object which contains the input data
+ including <type>BrookCorey</type>
+ and it has a tag of <capillary_pressure>
+*/
+static std::unique_ptr<CapillaryPressureSaturation> createBrookCorey(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__type}
+ config.checkConfigParameter("type", "BrookCorey");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__BrookCorey__pd}
+ const double pd = config.getConfigParameter<double>("pd");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__BrookCorey__sr}
+ const double Sr = config.getConfigParameter<double>("sr");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__BrookCorey__smax}
+ const double Smax = config.getConfigParameter<double>("smax");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__BrookCorey__m}
+ const double m = config.getConfigParameter<double>("m");
+ if (m < 1.0) // m >= 1
+ {
+ OGS_FATAL(
+ "The exponent parameter of BrookCorey capillary pressure "
+ "saturation model, m, must not be smaller than 1");
+ }
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__BrookCorey__pc_max}
+ const double Pc_max = config.getConfigParameter<double>("pc_max");
+
+ return std::unique_ptr<CapillaryPressureSaturation>(
+ new BrookCoreyCapillaryPressureSaturation(pd, Sr, Smax, m, Pc_max));
+}
+
+/**
+ \param config ConfigTree object which contains the input data
+ including <type>vanGenuchten</type>
+ and it has a tag of <capillary_pressure>
+*/
+static std::unique_ptr<CapillaryPressureSaturation> createVanGenuchten(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__type}
+ config.checkConfigParameter("type", "vanGenuchten");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__vanGenuchten__pd}
+ const double pd = config.getConfigParameter<double>("pd");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__vanGenuchten__sr}
+ const double Sr = config.getConfigParameter<double>("sr");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__vanGenuchten__smax}
+ const double Smax = config.getConfigParameter<double>("smax");
+
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__vanGenuchten__m}
+ const double m = config.getConfigParameter<double>("m");
+ if (m > 1.0) // m <= 1
+ {
+ OGS_FATAL(
+ "The exponent parameter of van Genuchten capillary pressure "
+ "saturation model, m, must be in an interval of [0, 1]");
+ }
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__vanGenuchten__pc_max}
+ const double Pc_max = config.getConfigParameter<double>("pc_max");
+
+ return std::unique_ptr<CapillaryPressureSaturation>(
+ new VanGenuchtenCapillaryPressureSaturation(pd, Sr, Smax, m, Pc_max));
+}
+
+std::unique_ptr<CapillaryPressureSaturation> createCapillaryPressureModel(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{material_property__porous_medium__porous_medium__capillary_pressure__type}
+ auto const type = config.peekConfigParameter<std::string>("type");
+
+ if (type == "BrookCorey")
+ {
+ return createBrookCorey(config);
+ }
+ else if (type == "vanGenuchten")
+ {
+ return createVanGenuchten(config);
+ }
+ else
+ {
+ OGS_FATAL(
+ "The capillary pressure models %s are unavailable.\n"
+ "The available types are: \n\tBrookCorey, \n\tvanGenuchten.\n",
+ type.data());
+ }
+}
+
+} // end namespace
+} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h
new file mode 100644
index 0000000000000000000000000000000000000000..fdedfc56d5b0688d5f68b4597377a59dceed2879
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h
@@ -0,0 +1,35 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file CreateCapillaryPressureModel.h
+ *
+ * Created on November 1, 2016, 10:06 AM
+ */
+
+#ifndef OGS_CREATE_CAPILLARY_PRESSURE_MODEL_H
+#define OGS_CREATE_CAPILLARY_PRESSURE_MODEL_H
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+class CapillaryPressureSaturation;
+/// Create a capillary pressure model
+/// \param config ConfigTree object has a tag of <capillary_pressure>
+std::unique_ptr<CapillaryPressureSaturation> createCapillaryPressureModel(
+ BaseLib::ConfigTree const& config);
+}
+} // end namespace
+
+#endif /* OGS_CREATE_CAPILLARY_PRESSURE_MODEL_H */
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c8cef0d9226b242fca3f94e58afd7dd8e369a346
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp
@@ -0,0 +1,57 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file VanGenuchtenCapillaryPressureSaturation.cpp
+ *
+ * Created on October 28, 2016, 6:05 PM
+ */
+
+#include "VanGenuchtenCapillaryPressureSaturation.h"
+
+#include <cmath>
+
+#include "MathLib/MathTools.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+double VanGenuchtenCapillaryPressureSaturation::getCapillaryPressure(
+ const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double Se = (S - _Sr) / (_Smax - _Sr);
+ const double pc =
+ _pb * std::pow(std::pow(Se, (-1.0 / _mm)) - 1.0, 1.0 - _mm);
+ return MathLib::limitValueInInterval(pc, _minor_offset, _Pc_max);
+}
+
+double VanGenuchtenCapillaryPressureSaturation::getSaturation(
+ const double capillary_pressure) const
+{
+ const double pc =
+ (capillary_pressure < 0.0) ? _minor_offset : capillary_pressure;
+ double Se = std::pow(pc / _pb, 1.0 / (1.0 - _mm)) + 1.0;
+ Se = std::pow(Se, -_mm);
+ const double S = Se * (_Smax - _Sr) + _Sr;
+ return MathLib::limitValueInInterval(S, _Sr + _minor_offset,
+ _Smax - _minor_offset);
+}
+
+double VanGenuchtenCapillaryPressureSaturation::getdPcdS(
+ const double saturation) const
+{
+ const double S = MathLib::limitValueInInterval(
+ saturation, _Sr + _minor_offset, _Smax - _minor_offset);
+ const double val1 = std::pow(((S - _Sr) / (_Smax - _Sr)), -1.0 / _mm);
+ const double val2 = std::pow(val1 - 1.0, -_mm);
+ return _pb * (_mm - 1.0) * val1 * val2 / (_mm * (S - _Sr));
+}
+
+} // end namespace
+} // end namespace
diff --git a/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h
new file mode 100644
index 0000000000000000000000000000000000000000..1b6fb3c33858b760efb32f684fa73486a13ca448
--- /dev/null
+++ b/MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h
@@ -0,0 +1,97 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file VanGenuchtenCapillaryPressureSaturation.h
+ *
+ * Created on October 28, 2016, 6:05 PM
+ */
+
+#ifndef OGS_VAN_GENUCHTEN_CAPILLARY_PRESSURE_SATURATION_H
+#define OGS_VAN_GENUCHTEN_CAPILLARY_PRESSURE_SATURATION_H
+
+#include <limits>
+
+#include "CapillaryPressureSaturation.h"
+
+namespace MaterialLib
+{
+namespace PorousMedium
+{
+/**
+ * \brief van Genuchten water retention model
+ *
+ * \f[p_c=p_b (S_e^{-1/m}-1)^{1-m}\f]
+ * with
+ * \f[S_e=\dfrac{S-S_r}{S_{\mbox{max}}-S_r}\f]
+ * where
+ * \f{eqnarray*}{
+ * &p_b& \mbox{ entry pressure,}\\
+ * &S_r& \mbox{ residual saturation,}\\
+ * &S_{\mbox{max}}& \mbox{ maximum saturation,}\\
+ * &m(<=1) & \mbox{ exponent.}\\
+ * \f}
+ *
+ * Note:
+ * \f[m=1/(1-n)\f].
+ *
+ * If \f$\alpha\f$ instead of \f$p_b\f$ is available, \f$p_b\f$ can be
+ * calculated
+ * as
+ * \f[p_b=\rho g/\alpha\f]
+ */
+class VanGenuchtenCapillaryPressureSaturation final
+ : public CapillaryPressureSaturation
+{
+public:
+ /**
+ * @param pb Entry pressure, \f$ p_b \f$
+ * @param Sr Residual saturation, \f$ S_r \f$
+ * @param Smax Maximum saturation, \f$ S_{\mbox{max}} \f$
+ * @param m Exponent, \f$ m \f$
+ * @param Pc_max Maximum capillary pressure, \f$ P_c^{\mbox{max}}\f$
+ */
+ VanGenuchtenCapillaryPressureSaturation(const double pb, const double Sr,
+ const double Smax, const double m,
+ const double Pc_max)
+ : _pb(pb), _Sr(Sr), _Smax(Smax), _mm(m), _Pc_max(Pc_max)
+ {
+ }
+
+ /// Get model name.
+ std::string getName() const override
+ {
+ return "van Genuchten water retention model.";
+ }
+
+ /// Get capillary pressure.
+ double getCapillaryPressure(const double saturation) const override;
+
+ /// Get capillary pressure.
+ double getSaturation(const double capillary_pressure) const override;
+
+ /// Get the derivative of the capillary pressure with respect to saturation
+ double getdPcdS(const double saturation) const override;
+
+private:
+ const double _pb; ///< Entry pressure.
+ const double _Sr; ///< Residual saturation.
+ const double _Smax; ///< Maximum saturation.
+ const double _mm; ///< Exponent (<=1.0), n=1/(1-mm).
+ const double _Pc_max; ///< Maximum capillaray pressure
+
+ /** A small number for an offset:
+ * 1. to set the bound of S, the saturation, such that
+ * S in [_Sr+_minor_offset, _Smax-_minor_offset]
+ * 2. to set the bound of Pc, the capillary pressure, such that
+ * Pc in [_minor_offset, _Pc_max]
+ */
+ const double _minor_offset = std::numeric_limits<double>::epsilon();
+};
+
+} // end namespace
+} // end namespace
+#endif /* OGS_VAN_GENUCHTEN_CAPILLARY_PRESSURE_SATURATION_H */
diff --git a/MathLib/MathTools.h b/MathLib/MathTools.h
index e578be54451e0447bbd9aa198bf93c8c2c21140e..5f37210fffa8593d00b5032ff9af43399bdabb25 100644
--- a/MathLib/MathTools.h
+++ b/MathLib/MathTools.h
@@ -120,6 +120,16 @@ inline double to_radians(double degrees) {
return degrees*boost::math::constants::pi<double>()/180.;
}
+template<typename Type> Type limitValueInInterval(const Type variable,
+ const Type lower_bound,
+ const Type upper_bound)
+{
+ if (variable < lower_bound)
+ return lower_bound;
+ if (variable > upper_bound)
+ return upper_bound;
+ return variable;
+}
} // namespace
#endif /* MATHTOOLS_H_ */
diff --git a/Tests/MaterialLib/TestCapillaryPressureSaturationModel.cpp b/Tests/MaterialLib/TestCapillaryPressureSaturationModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..98d9d7c1fe99e2625982085888e5579cc6221acc
--- /dev/null
+++ b/Tests/MaterialLib/TestCapillaryPressureSaturationModel.cpp
@@ -0,0 +1,103 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * \file TestCapillaryPressureSaturationModel.cpp
+ *
+ * Created on November 1, 2016, 11:06 AM
+ */
+
+#include <gtest/gtest.h>
+
+#include <memory>
+#include <vector>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "TestTools.h"
+
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h"
+#include "MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h"
+
+using namespace MaterialLib;
+using namespace MaterialLib::PorousMedium;
+
+std::unique_ptr<CapillaryPressureSaturation> createCapillaryPressureModel(
+ const char xml[])
+{
+ auto const ptree = readXml(xml);
+ BaseLib::ConfigTree conf(ptree, "", BaseLib::ConfigTree::onerror,
+ BaseLib::ConfigTree::onwarning);
+ auto const& sub_config = conf.getConfigSubtree("capillary_pressure");
+ return MaterialLib::PorousMedium::createCapillaryPressureModel(sub_config);
+}
+
+TEST(MaterialPorousMedium, checkBrookCoreyCapillaryPressure)
+{
+ const char xml[] =
+ "<capillary_pressure>"
+ " <type>BrookCorey</type>"
+ " <pd> 19600.0 </pd> "
+ " <sr> 0.1 </sr> "
+ " <smax> 1. </smax> "
+ " <m> 2 </m> "
+ " <pc_max> 1.e11 </pc_max> "
+ "</capillary_pressure>";
+ auto const pc_model = createCapillaryPressureModel(xml);
+
+ std::vector<double> S = {0.11, 0.2, 0.3, 0.44, 0.52, 0.6, 0.85};
+ // The following expected data are calculated by using OGS5.
+ std::vector<double> pc = {185941.926417901, 58800,
+ 41577.878733769, 31888.7772993424,
+ 28691.4621446869, 26296.1594153975,
+ 21470.7242542025};
+ std::vector<double> dpc_dS = {-9297096.32089506, -294000,
+ -103944.696834422, -46895.2607343271,
+ -34156.5025531986, -26296.1594153975,
+ -14313.8161694683};
+
+ const double tol_pc = 1.e-6;
+ for (std::size_t i = 0; i < S.size(); i++)
+ {
+ ASSERT_NEAR(S[i], pc_model->getSaturation(pc[i]), 1e-14);
+ ASSERT_NEAR(pc[i], pc_model->getCapillaryPressure(S[i]), tol_pc);
+ ASSERT_NEAR(dpc_dS[i], pc_model->getdPcdS(S[i]), tol_pc);
+ }
+}
+
+// In the following test, the expected results are calculated by using OGS5.
+TEST(MaterialPorousMedium, checkVanGenuchtenCapillaryPressure)
+{
+ // rho=1000, alpha = 0.37, pd=rho*g/alpha
+ const char xml[] =
+ "<capillary_pressure>"
+ " <type>vanGenuchten</type>"
+ " <pd> 26513.513513513513 </pd> "
+ " <sr> 0. </sr> "
+ " <smax> 1. </smax> "
+ " <m> 0.7 </m> "
+ " <pc_max> 1.e5 </pc_max> "
+ "</capillary_pressure>";
+ auto const pc_model = createCapillaryPressureModel(xml);
+
+ std::vector<double> S = {0.11, 0.2, 0.3, 0.44, 0.52, 0.6, 0.85};
+ // The following expected data are calculated by using OGS5.
+ std::vector<double> pc = {
+ 67392.5090996789, 51197.5842770154, 41864.8480636163, 33730.6152475992,
+ 30210.4060771713, 27091.7425032625, 17726.625235802};
+ std::vector<double> dpc_dS = {-274283.722262232, -121944.994573743,
+ -72852.9253517949, -47580.0457562877,
+ -41013.0131889419, -37359.7105247455,
+ -43138.4488851645};
+
+ const double tol_pc = 1.e-6;
+ for (std::size_t i = 0; i < S.size(); i++)
+ {
+ ASSERT_NEAR(S[i], pc_model->getSaturation(pc[i]), 1e-14);
+ ASSERT_NEAR(pc[i], pc_model->getCapillaryPressure(S[i]), tol_pc);
+ ASSERT_NEAR(dpc_dS[i], pc_model->getdPcdS(S[i]), tol_pc);
+ }
+}