Merge pull request #1517 from wenqing/perm_Pc_S

Capillary pressure models

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

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.cpp

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

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/BrookCoreyCapillaryPressureSaturation.h

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

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CapillaryPressureSaturation.h

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

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.cpp

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

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/CreateCapillaryPressureModel.h

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

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.cpp

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

MaterialLib/PorousMedium/UnsaturatedProperty/CapillaryPressure/VanGenuchtenCapillaryPressureSaturation.h

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

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_ */

Tests/MaterialLib/TestCapillaryPressureSaturationModel.cpp

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