Commit 76f0149b authored by Dmitry Yu. Naumov's avatar Dmitry Yu. Naumov
Browse files

Merge branch 'MPL_SaturationExponential' into 'master'

MPL SaturationExponential

See merge request ogs/ogs!3579
parents 1594b566 001cb341
\copydoc MaterialPropertyLib::SaturationExponential::exponent_
\copydoc MaterialPropertyLib::SaturationExponential::residual_gas_saturation_
\copydoc MaterialPropertyLib::SaturationExponential::residual_liquid_saturation_
......@@ -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);
......
/**
* \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_max =
//! \ogs_file_param{properties__property__SaturationExponential__maximum_capillary_pressure}
config.getConfigParameter<double>("maximum_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_max, exponent);
}
} // namespace MaterialPropertyLib
/**
* \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
/**
* \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_max,
const double exponent)
: residual_liquid_saturation_(residual_liquid_saturation),
residual_gas_saturation_(residual_gas_saturation),
p_cap_max_(p_cap_max),
exponent_(exponent)
{
name_ = std::move(name);
if (p_cap_max_ <= 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_max_);
const double S_e = 1. - std::pow(pc / p_cap_max_, 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_max_) || (p_cap <= 0.))
{
return 0.;
}
return (exponent_ / p_cap) * (S_res - S_max) *
std::pow(p_cap / p_cap_max_, 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
/**
* \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 "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.
*
* The capillary pressure-saturation relationship given by
* \f$S_{eff}=1-\left(\frac{p_{c}}{p_{c}^{max}}\right)^{\lambda}\f$,
* where
* \f$\lambda\f$ is an exponent,
* \f$p_{c}\f$ is capillary pressure,
* \f$p_{c}^{max}\f$ is the maximum capillary pressure at which \f$S_{eff}=0\f$.
*
* This property can mainly be used for testing. If the exponent is set to 1,
* the characteristic curve shows a linear dependence.
*
*
* Reference capillary pressure at which \f$S_{eff}=0\f$.
*/
class SaturationExponential final : public Property
{
private:
const double
residual_liquid_saturation_; ///< Residual saturation of the gas phase.
const double
residual_gas_saturation_; ///< Residual saturation of the liquid phase.
const double p_cap_max_; ///< Maximum capillary pressure at which
/// effective saturation is zero.
const double exponent_; ///< Exponent to govern the shape of the curve.
public:
SaturationExponential(std::string name,
const double residual_liquid_saturation,
const double residual_gas_saturation,
const double p_cap_max, 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.");
}
}
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
......@@ -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"
......
......@@ -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"
......
/**
* \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 <gtest/gtest.h>
#include <iomanip>
#include "MaterialLib/MPL/Medium.h"
#include "MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h"
#include "TestMPL.h"
#include "Tests/TestTools.h"
namespace MPL = MaterialPropertyLib;
TEST(MaterialPropertyLib, SaturationExponential)
{
double const residual_liquid_saturation = 0.1;
double const residual_gas_saturation = 0.05;
double const exponent = 1.5;
double const p_cap_max = 30000;
MPL::Property const& saturation = MPL::SaturationExponential{
"saturation", residual_liquid_saturation, residual_gas_saturation,
p_cap_max, exponent};
MPL::VariableArray variable_array;
ParameterLib::SpatialPosition const pos;
double const t = std::numeric_limits<double>::quiet_NaN();
double const dt = std::numeric_limits<double>::quiet_NaN();
double const p_0 = -p_cap_max;
double const p_max = 1.5 * p_cap_max;
int const n_steps = 9999;
for (int i = 0; i <= n_steps; ++i)
{
double const p_cap = p_0 + i * (p_max - p_0) / n_steps;
variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap;
double const s_res = residual_liquid_saturation;
double const s_max = 1. - residual_gas_saturation;
double const s_eff =
1. -
std::pow(std::clamp(p_cap, 0., p_cap_max) / p_cap_max, exponent);
double s_ref = s_eff * (s_max - s_res) + s_res;
double const S =
saturation.template value<double>(variable_array, pos, t, dt);
double const dS = saturation.template dValue<double>(
variable_array, MPL::Variable::capillary_pressure, pos, t, dt);
double const eps = 1e-2;
variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap - eps;
double const S_minus =
saturation.template value<double>(variable_array, pos, t, dt);
variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap + eps;
double const S_plus =
saturation.template value<double>(variable_array, pos, t, dt);
double const dS_ref = p_cap > 0 ? (S_plus - S_minus) / 2 / eps : 0.;
ASSERT_LE(std::abs(s_ref - S), 1e-9)
<< std::setprecision(16) << "with: \ncapillary pressure: " << p_cap
<< "\nsaturation: " << S << "\nsaturation_ref = " << s_ref
<< "\nat point: " << i;
ASSERT_LE(std::abs(dS - dS_ref), 1e-9)
<< std::setprecision(16) << "with: \np_cap: " << p_cap
<< "\nS: " << S << "\nS_ref = " << s_ref << "\ndS = " << dS
<< "\ndS_ref = " << dS_ref << "\nat point: " << i;
}
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment