<MPL> add saturationExponential property

MaterialLib/MPL/CreateProperty.cpp

@@ -157,6 +157,10 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
         return createRelPermLiakopoulos(config);
     }
 
+    if (boost::iequals(property_type, "SaturationExponential"))
+    {
+        return createSaturationExponential(config);
+    }
     if (boost::iequals(property_type, "SaturationVanGenuchten"))
     {
         return createSaturationVanGenuchten(config);

MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "SaturationExponential.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationExponential> createSaturationExponential(
+    BaseLib::ConfigTree const& config)
+{
+    //! \ogs_file_param{properties__property__type}
+    config.checkConfigParameter("type", "SaturationExponential");
+
+    // Second access for storage.
+    //! \ogs_file_param{properties__property__name}
+    auto property_name = config.peekConfigParameter<std::string>("name");
+
+    DBUG("Create SaturationExponential medium property {:s}.", property_name);
+
+    auto const residual_liquid_saturation =
+        //! \ogs_file_param{properties__property__SaturationExponential__residual_liquid_saturation}
+        config.getConfigParameter<double>("residual_liquid_saturation");
+    auto const residual_gas_saturation =
+        //! \ogs_file_param{properties__property__SaturationExponential__residual_gas_saturation}
+        config.getConfigParameter<double>("residual_gas_saturation");
+    auto const p_cap_ref =
+        //! \ogs_file_param{properties__property__SaturationExponential__reference_capillary_pressure}
+        config.getConfigParameter<double>("reference_capillary_pressure");
+    auto const exponent =
+        //! \ogs_file_param{properties__property__SaturationExponential__exponent}
+        config.getConfigParameter<double>("exponent");
+
+    return std::make_unique<SaturationExponential>(
+        std::move(property_name), residual_liquid_saturation,
+        residual_gas_saturation, p_cap_ref, exponent);
+}
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h

+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class SaturationExponential;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationExponential> createSaturationExponential(
+    BaseLib::ConfigTree const& config);
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp

+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date   27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "SaturationExponential.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MathLib/MathTools.h"
+
+namespace MaterialPropertyLib
+{
+SaturationExponential::SaturationExponential(
+    std::string name,
+    const double residual_liquid_saturation,
+    const double residual_gas_saturation,
+    const double p_cap_ref,
+    const double exponent)
+    : residual_liquid_saturation_(residual_liquid_saturation),
+      residual_gas_saturation_(residual_gas_saturation),
+      p_cap_ref_(p_cap_ref),
+      exponent_(exponent)
+{
+    name_ = std::move(name);
+    if (p_cap_ref_ <= 0.)
+    {
+        OGS_FATAL(
+            "Reference capillary pressure must be positive in "
+            "MPL::SaturationExponential.");
+    }
+};
+
+PropertyDataType SaturationExponential::value(
+    VariableArray const& variable_array,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    const double p_cap = std::get<double>(
+        variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+    const double s_res = residual_liquid_saturation_;
+    const double s_max = 1. - residual_gas_saturation_;
+
+    const double pc = std::clamp(p_cap, 0., p_cap_ref_);
+    const double s_e = 1. - std::pow(pc / p_cap_ref_, exponent_);
+    return s_e * (s_max - s_res) + s_res;
+}
+
+PropertyDataType SaturationExponential::dValue(
+    VariableArray const& variable_array, Variable const primary_variable,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    (void)primary_variable;
+    assert((primary_variable == Variable::capillary_pressure) &&
+           "SaturationExponential::dValue is implemented for derivatives with "
+           "respect to capillary pressure only.");
+
+    const double p_cap = std::get<double>(
+        variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+    const double s_res = residual_liquid_saturation_;
+    const double s_max = 1. - residual_gas_saturation_;
+
+    if ((p_cap > p_cap_ref_) || (p_cap <= 0.))
+    {
+        return 0.;
+    }
+    return (exponent_ / p_cap) * (s_res - s_max) *
+           std::pow(p_cap / p_cap_ref_, exponent_);
+}
+
+PropertyDataType SaturationExponential::d2Value(
+    VariableArray const& /*variable_array*/,
+    Variable const /*primary_variable1*/, Variable const /*primary_variable2*/,
+    ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+    double const /*dt*/) const
+{
+    OGS_FATAL("SaturationExponential::d2Value() is not implemented.");
+}
+
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h

+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date   27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, 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
+{
+/**
+ * \brief A simplistic 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 SaturationExponential final : public Property
+{
+private:
+    const double residual_liquid_saturation_;
+    const double residual_gas_saturation_;
+    const double p_cap_ref_;
+    const double exponent_;
+
+public:
+    SaturationExponential(std::string name,
+                          const double residual_liquid_saturation,
+                          const double residual_gas_saturation,
+                          const double p_cap_ref, const double exponent);
+
+    void checkScale() const override
+    {
+        if (!std::holds_alternative<Medium*>(scale_))
+        {
+            OGS_FATAL(
+                "The property 'SaturationExponential' is implemented on the "
+                "'media' scale only.");
+        }
+    }
+
+    /// 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*/,
+                           double const /*dt*/) const override;
+    PropertyDataType dValue(VariableArray const& variable_array,
+                            Variable const primary_variable,
+                            ParameterLib::SpatialPosition const& /*pos*/,
+                            double const /*t*/,
+                            double const /*dt*/) const override;
+    PropertyDataType d2Value(VariableArray const& variable_array,
+                             Variable const /*primary_variable1*/,
+                             Variable const /*primary_variable2*/,
+                             ParameterLib::SpatialPosition const& /*pos*/,
+                             double const /*t*/,
+                             double const /*dt*/) const override;
+};
+}  // namespace MaterialPropertyLib

MaterialLib/MPL/Properties/CreateProperties.h

@@ -15,6 +15,7 @@
 #include "CapillaryPressureSaturation/CreateCapillaryPressureVanGenuchten.h"
 #include "CapillaryPressureSaturation/CreateSaturationBrooksCorey.h"
 #include "CapillaryPressureSaturation/CreateSaturationLiakopoulos.h"
+#include "CapillaryPressureSaturation/CreateSaturationExponential.h"
 #include "CapillaryPressureSaturation/CreateSaturationVanGenuchten.h"
 #include "CreateAverageMolarMass.h"
 #include "CreateBishopsPowerLaw.h"

MaterialLib/MPL/Properties/Properties.h

@@ -16,6 +16,7 @@
 #include "BishopsSaturationCutoff.h"
 #include "CapillaryPressureSaturation/SaturationBrooksCorey.h"
 #include "CapillaryPressureSaturation/SaturationLiakopoulos.h"
+#include "CapillaryPressureSaturation/SaturationExponential.h"
 #include "CapillaryPressureSaturation/SaturationVanGenuchten.h"
 #include "ClausiusClapeyron.h"
 #include "Constant.h"