From 8a455bb79b09f9e71ddfa87e8da5a5e1a27c7260 Mon Sep 17 00:00:00 2001
From: Norbert Grunwald <Norbert.Grunwald@ufz.de>
Date: Mon, 9 Nov 2020 14:13:13 +0100
Subject: [PATCH] <MPL> add saturationExponential property
---
MaterialLib/MPL/CreateProperty.cpp | 4 +
.../CreateSaturationExponential.cpp | 45 +++++++++
.../CreateSaturationExponential.h | 29 ++++++
.../SaturationExponential.cpp | 94 +++++++++++++++++++
.../SaturationExponential.h | 70 ++++++++++++++
MaterialLib/MPL/Properties/CreateProperties.h | 1 +
MaterialLib/MPL/Properties/Properties.h | 1 +
7 files changed, 244 insertions(+)
create mode 100644 MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp
create mode 100644 MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h
create mode 100644 MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp
create mode 100644 MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 049f4a7e1a2..082de71d9be 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/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);
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp
new file mode 100644
index 00000000000..f0778c81c38
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp
@@ -0,0 +1,45 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h
new file mode 100644
index 00000000000..2e00abdee2f
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h
@@ -0,0 +1,29 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp
new file mode 100644
index 00000000000..9e35916b4b1
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp
@@ -0,0 +1,94 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h
new file mode 100644
index 00000000000..3f1063d687f
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h
@@ -0,0 +1,70 @@
+/**
+ * \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
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 7f177a085f4..4c62b2242dd 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/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"
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index ee5a30636d4..f1ec5467d48 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/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"
--
GitLab