Loading...
Searching...
No Matches
EclDefaultMaterial.hpp
Go to the documentation of this file.
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
#ifndef OPM_ECL_DEFAULT_MATERIAL_HPP
#define OPM_ECL_DEFAULT_MATERIAL_HPP
#include <opm/common/TimingMacros.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp>
#include <algorithm>
#include <stdexcept>
#include <type_traits>
namespace Opm {
template <class TraitsT,
class GasOilMaterialLawT,
class OilWaterMaterialLawT,
class ParamsT = EclDefaultMaterialParams<TraitsT,
typename GasOilMaterialLawT::Params,
typename OilWaterMaterialLawT::Params> >
class EclDefaultMaterial : public TraitsT
{
public:
using GasOilMaterialLaw = GasOilMaterialLawT;
using OilWaterMaterialLaw = OilWaterMaterialLawT;
// some safety checks
static_assert(TraitsT::numPhases == 3,
"The number of phases considered by this capillary pressure "
"law is always three!");
static_assert(GasOilMaterialLaw::numPhases == 2,
"The number of phases considered by the gas-oil capillary "
"pressure law must be two!");
static_assert(OilWaterMaterialLaw::numPhases == 2,
"The number of phases considered by the oil-water capillary "
"pressure law must be two!");
static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
typename OilWaterMaterialLaw::Scalar>::value,
"The two two-phase capillary pressure laws must use the same "
"type of floating point values.");
static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
"The gas-oil material law must implement the two-phase saturation "
"only API to for the default Ecl capillary pressure law!");
static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
"The oil-water material law must implement the two-phase saturation "
"only API to for the default Ecl capillary pressure law!");
using Traits = TraitsT;
using Params = ParamsT;
using Scalar = typename Traits::Scalar;
static constexpr int numPhases = 3;
static constexpr int waterPhaseIdx = Traits::wettingPhaseIdx;
static constexpr int oilPhaseIdx = Traits::nonWettingPhaseIdx;
static constexpr int gasPhaseIdx = Traits::gasPhaseIdx;
static constexpr bool implementsTwoPhaseApi = false;
static constexpr bool implementsTwoPhaseSatApi = false;
static constexpr bool isSaturationDependent = true;
static constexpr bool isPressureDependent = false;
static constexpr bool isTemperatureDependent = false;
static constexpr bool isCompositionDependent = false;
template <class ContainerT, class FluidState, class ...Args>
static void capillaryPressures(ContainerT& values,
const Params& params,
const FluidState& state)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
values[gasPhaseIdx] = pcgn<FluidState, Evaluation, Args...>(params, state);
values[oilPhaseIdx] = 0;
values[waterPhaseIdx] = - pcnw<FluidState, Evaluation, Args...>(params, state);
Valgrind::CheckDefined(values[gasPhaseIdx]);
Valgrind::CheckDefined(values[oilPhaseIdx]);
Valgrind::CheckDefined(values[waterPhaseIdx]);
}
/*
* Hysteresis parameters for oil-water
* @see EclHysteresisTwoPhaseLawParams::soMax(...)
* @see EclHysteresisTwoPhaseLawParams::swMax(...)
* @see EclHysteresisTwoPhaseLawParams::swMin(...)
* \param params Parameters
*/
static void oilWaterHysteresisParams(Scalar& soMax,
Scalar& swMax,
Scalar& swMin,
const Params& params)
{
if constexpr (Traits::enableHysteresis) {
soMax = 1.0 - params.oilWaterParams().krnSwMdc();
swMax = params.oilWaterParams().krwSwMdc();
swMin = params.oilWaterParams().pcSwMdc();
Valgrind::CheckDefined(soMax);
Valgrind::CheckDefined(swMax);
Valgrind::CheckDefined(swMin);
}
}
/*
* Hysteresis parameters for oil-water
* @see EclHysteresisTwoPhaseLawParams::soMax(...)
* @see EclHysteresisTwoPhaseLawParams::swMax(...)
* @see EclHysteresisTwoPhaseLawParams::swMin(...)
* \param params Parameters
*/
static void setOilWaterHysteresisParams(const Scalar& soMax,
const Scalar& swMax,
const Scalar& swMin,
Params& params)
{
if constexpr (Traits::enableHysteresis) {
params.oilWaterParams().update(swMin, swMax, 1.0 - soMax);
}
}
/*
* Hysteresis parameters for gas-oil
* @see EclHysteresisTwoPhaseLawParams::sgMax(...)
* @see EclHysteresisTwoPhaseLawParams::shMax(...)
* @see EclHysteresisTwoPhaseLawParams::soMin(...)
* \param params Parameters
*/
static void gasOilHysteresisParams(Scalar& sgMax,
Scalar& shMax,
Scalar& soMin,
const Params& params)
{
if constexpr (Traits::enableHysteresis) {
const auto Swco = params.Swl();
sgMax = 1.0 - params.gasOilParams().krnSwMdc() - Swco;
shMax = params.gasOilParams().krwSwMdc();
soMin = params.gasOilParams().pcSwMdc();
Valgrind::CheckDefined(sgMax);
Valgrind::CheckDefined(shMax);
Valgrind::CheckDefined(soMin);
}
}
/*
* Hysteresis parameters for gas-oil
* @see EclHysteresisTwoPhaseLawParams::sgMax(...)
* @see EclHysteresisTwoPhaseLawParams::shMax(...)
* @see EclHysteresisTwoPhaseLawParams::soMin(...)
* \param params Parameters
*/
static void setGasOilHysteresisParams(const Scalar& sgMax,
const Scalar& shMax,
const Scalar& soMin,
Params& params)
{
if constexpr (Traits::enableHysteresis) {
const auto Swco = params.Swl();
params.gasOilParams().update(soMin, shMax, 1.0 - sgMax - Swco);
}
}
static Scalar trappedGasSaturation(const Params& params, bool maximumTrapping)
{
const auto Swco = params.Swl();
return params.gasOilParams().SnTrapped(maximumTrapping) - Swco;
}
static Scalar trappedOilSaturation(const Params& params, bool maximumTrapping)
{
return params.oilWaterParams().SnTrapped(maximumTrapping) + params.gasOilParams().SwTrapped();
}
static Scalar trappedWaterSaturation(const Params& params)
{
return params.oilWaterParams().SwTrapped();
}
static Scalar strandedGasSaturation(const Params& params, Scalar Sg, Scalar Kg)
{
const auto Swco = params.Swl();
return params.gasOilParams().SnStranded(Sg, Kg) - Swco;
}
template <class FluidState, class Evaluation, class ...Args>
static Evaluation pcgn(const Params& params,
const FluidState& fs)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
// Maximum attainable oil saturation is 1-SWL.
const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::template twoPhaseSatPcnw<Evaluation, Args...>(params.gasOilParams(), Sw);
}
template <class FluidState, class Evaluation, class ...Args>
static Evaluation pcnw(const Params& params,
const FluidState& fs)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::template twoPhaseSatPcnw<Evaluation, Args...>(params.oilWaterParams(), Sw);
}
template <class ContainerT, class FluidState>
static void saturations(ContainerT& /*values*/,
const Params& /*params*/,
const FluidState& /*fluidState*/)
{
throw std::logic_error("Not implemented: saturations()");
}
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sg(const Params& /*params*/,
const FluidState& /*fluidState*/)
{
throw std::logic_error("Not implemented: Sg()");
}
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params& /*params*/,
const FluidState& /*fluidState*/)
{
throw std::logic_error("Not implemented: Sn()");
}
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params& /*params*/,
const FluidState& /*fluidState*/)
{
throw std::logic_error("Not implemented: Sw()");
}
template <class ContainerT, class FluidState, class ...Args>
static void relativePermeabilities(ContainerT& values,
const Params& params,
const FluidState& fluidState)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
values[waterPhaseIdx] = krw<FluidState, Evaluation, Args...>(params, fluidState);
values[oilPhaseIdx] = krn<FluidState, Evaluation, Args...>(params, fluidState);
values[gasPhaseIdx] = krg<FluidState, Evaluation, Args...>(params, fluidState);
}
template <class FluidState, class Evaluation, class ...Args>
static Evaluation krg(const Params& params,
const FluidState& fluidState)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
// Maximum attainable oil saturation is 1-SWL.
const Evaluation sw = 1.0 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::template twoPhaseSatKrn<Evaluation, Args...>(params.gasOilParams(), sw);
}
template <class FluidState, class Evaluation, class ...Args>
static Evaluation krw(const Params& params,
const FluidState& fluidState)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
const Evaluation sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::template twoPhaseSatKrw<Evaluation, Args...>(params.oilWaterParams(), sw);
}
template <class FluidState, class Evaluation, class ...Args>
static Evaluation krn(const Params& params,
const FluidState& fluidState)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
const Scalar Swco = params.Swl();
const Evaluation sw =
max(Evaluation(Swco),
decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
const Evaluation sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation Sw_ow = sg + sw;
const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation, FluidState, Args...>(params, fluidState);
const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation, FluidState, Args...>(params, fluidState);
// avoid the division by zero: chose a regularized kro which is used if Sw - Swco
// < epsilon/2 and interpolate between the oridinary and the regularized kro between
// epsilon and epsilon/2
constexpr const Scalar epsilon = 1e-5;
if (scalarValue(Sw_ow) - Swco < epsilon) {
const Evaluation kro2 = (kro_ow + kro_go)/2;
if (scalarValue(Sw_ow) - Swco > epsilon/2) {
const Evaluation kro1 = (sg * kro_go + (sw - Swco) * kro_ow) / (Sw_ow - Swco);
const Evaluation alpha = (epsilon - (Sw_ow - Swco)) / (epsilon / 2);
return kro2 * alpha + kro1 * (1 - alpha);
}
return kro2;
}
return (sg * kro_go + (sw - Swco) * kro_ow) / (Sw_ow - Swco);
}
template <class Evaluation, class FluidState, class ...Args>
static Evaluation relpermOilInOilGasSystem(const Params& params,
const FluidState& fluidState)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
const Evaluation sw =
max(Evaluation{ params.Swl() },
decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
const Evaluation sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation So_go = 1.0 - (sg + sw);
return GasOilMaterialLaw::template twoPhaseSatKrw<Evaluation, Args...>(params.gasOilParams(), So_go);
}
template <class Evaluation, class FluidState, class ...Args>
static Evaluation relpermOilInOilWaterSystem(const Params& params,
const FluidState& fluidState)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
const Evaluation sw =
max(Evaluation{ params.Swl() },
decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
const Evaluation sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation Sw_ow = sg + sw;
return OilWaterMaterialLaw::template twoPhaseSatKrn<Evaluation, Args...>(params.oilWaterParams(), Sw_ow);
}
template <class FluidState>
static bool updateHysteresis(Params& params, const FluidState& fluidState)
{
if constexpr (Traits::enableHysteresis) {
const Scalar Swco = params.Swl();
const Scalar sw = clampSaturation(fluidState, waterPhaseIdx);
const Scalar So = clampSaturation(fluidState, oilPhaseIdx);
const Scalar sg = clampSaturation(fluidState, gasPhaseIdx);
bool owChanged = params.oilWaterParams().update(/*pcSw=*/sw, /*krwSw=*/sw, /*krnSw=*/1 - So);
bool gochanged = params.gasOilParams().update(/*pcSw=*/So,
/*krwSw=*/So,
/*krnSw=*/1.0 - Swco - sg);
return owChanged || gochanged;
} else {
return false;
}
}
template <class FluidState>
static Scalar clampSaturation(const FluidState& fluidState, const int phaseIndex)
{
OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
const auto sat = scalarValue(fluidState.saturation(phaseIndex));
return std::clamp(sat, Scalar{0.0}, Scalar{1.0});
}
};
} // namespace Opm
#endif
EclDefaultMaterialParams.hpp
Default implementation for the parameters required by the default three-phase capillary pressure mode...
MathToolbox.hpp
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
Valgrind.hpp
Some templates to wrap the valgrind client request macros.
Opm::Valgrind::CheckDefined
OPM_HOST_DEVICE bool CheckDefined(const T &value)
Make valgrind complain if any of the memory occupied by an object is undefined.
Definition Valgrind.hpp:76
Opm::EclDefaultMaterialParams
Default implementation for the parameters required by the default three-phase capillary pressure mode...
Definition EclDefaultMaterialParams.hpp:46
Opm::EclDefaultMaterialParams::Swl
Scalar Swl() const
Return the saturation of "connate" water.
Definition EclDefaultMaterialParams.hpp:115
Opm::EclDefaultMaterialParams::oilWaterParams
const OilWaterParams & oilWaterParams() const
The parameter object for the oil-water twophase law.
Definition EclDefaultMaterialParams.hpp:83
Opm::EclDefaultMaterialParams::gasOilParams
const GasOilParams & gasOilParams() const
The parameter object for the gas-oil twophase law.
Definition EclDefaultMaterialParams.hpp:65
Opm::EclDefaultMaterial
Implements the default three phase capillary pressure law used by the ECLipse simulator.
Definition EclDefaultMaterial.hpp:62
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::krn
static Evaluation krn(const Params &params, const FluidState &fluidState)
Definition EclDefaultMaterial.hpp:389
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::relpermOilInOilWaterSystem
static Evaluation relpermOilInOilWaterSystem(const Params &params, const FluidState &fluidState)
Definition EclDefaultMaterial.hpp:446
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::isCompositionDependent
static constexpr bool isCompositionDependent
Definition EclDefaultMaterial.hpp:120
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::pcnw
static Evaluation pcnw(const Params &params, const FluidState &fs)
Definition EclDefaultMaterial.hpp:283
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::isTemperatureDependent
static constexpr bool isTemperatureDependent
Definition EclDefaultMaterial.hpp:116
Opm::EclDefaultMaterial::Sg
static Evaluation Sg(const Params &, const FluidState &)
The saturation of the gas phase.
Definition EclDefaultMaterial.hpp:306
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::implementsTwoPhaseApi
static constexpr bool implementsTwoPhaseApi
Definition EclDefaultMaterial.hpp:100
Opm::EclDefaultMaterial::saturations
static void saturations(ContainerT &, const Params &, const FluidState &)
The inverse of the capillary pressure.
Definition EclDefaultMaterial.hpp:295
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::relpermOilInOilGasSystem
static Evaluation relpermOilInOilGasSystem(const Params &params, const FluidState &fluidState)
Definition EclDefaultMaterial.hpp:428
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::krw
static Evaluation krw(const Params &params, const FluidState &fluidState)
Definition EclDefaultMaterial.hpp:377
Opm::EclDefaultMaterial::Sn
static Evaluation Sn(const Params &, const FluidState &)
The saturation of the non-wetting (i.e., oil) phase.
Definition EclDefaultMaterial.hpp:316
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::isSaturationDependent
static constexpr bool isSaturationDependent
Definition EclDefaultMaterial.hpp:108
Opm::EclDefaultMaterial::updateHysteresis
static bool updateHysteresis(Params &params, const FluidState &fluidState)
Update the hysteresis parameters after a time step.
Definition EclDefaultMaterial.hpp:468
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::krg
static Evaluation krg(const Params &params, const FluidState &fluidState)
Definition EclDefaultMaterial.hpp:364
Opm::EclDefaultMaterial::capillaryPressures
static void capillaryPressures(ContainerT &values, const Params &params, const FluidState &state)
Implements the default three phase capillary pressure law used by the ECLipse simulator.
Definition EclDefaultMaterial.hpp:137
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::pcgn
static Evaluation pcgn(const Params &params, const FluidState &fs)
Definition EclDefaultMaterial.hpp:264
Opm::EclDefaultMaterial::relativePermeabilities
static void relativePermeabilities(ContainerT &values, const Params &params, const FluidState &fluidState)
The relative permeability of all phases.
Definition EclDefaultMaterial.hpp:348
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::Sw
static Evaluation Sw(const Params &, const FluidState &)
Definition EclDefaultMaterial.hpp:326
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::isPressureDependent
static constexpr bool isPressureDependent
Definition EclDefaultMaterial.hpp:112
Opm::EclDefaultMaterial< TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw >::implementsTwoPhaseSatApi
static constexpr bool implementsTwoPhaseSatApi
Definition EclDefaultMaterial.hpp:104
Opm
This class implements a small container which holds the transmissibility mulitpliers for all the face...
Definition Exceptions.hpp:30