OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

make the fluid-matrix interactions usable with the local-AD framework

Browse Source 
also, adapt the unit test to test this capability.
This commit is contained in:
Andreas Lauser 2015-05-21 15:33:18 +02:00
parent 26c8553902
commit 8aaf2bfdab
17 changed files with 1124 additions and 1006 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

186
opm/material/fluidmatrixinteractions/BrooksCorey.hpp
View File

@ -27,6 +27,8 @@
#include "BrooksCoreyParams.hpp"
#include <opm/material/common/MathToolbox.hpp>
#include <algorithm>
#include <cassert>
#include <cmath>
@ -92,8 +94,10 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
}
/*!
@ -103,7 +107,9 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = Sw(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
@ -120,8 +126,10 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
}
/*!
@ -137,18 +145,27 @@ public:
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
assert(0 <= Sw && Sw <= 1);
return twoPhaseSatPcnw(params, Sw);
}
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0 <= Sw && Sw <= 1);
return params.entryPressure()*std::pow(Sw, -1.0/params.lambda());
return params.entryPressure()*Toolbox::pow(Sw, -1.0/params.lambda());
}
/*!
@ -163,78 +180,38 @@ public:
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{
Scalar pc = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
return twoPhaseSatSw(params, pc);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Evaluation pC =
FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pc)
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pc)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(pc > 0); // if we don't assume that, std::pow will screw up!
return std::pow(pc/params.entryPressure(), -params.lambda());
return Toolbox::pow(pc/params.entryPressure(), -params.lambda());
}
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw<FluidState, Evaluation>(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pc)
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pc)
{ return 1 - twoPhaseSatSw(params, pc); }
/*!
* \brief The partial derivative of the capillary
* pressure w.r.t. the effective saturation according to Brooks & Corey.
*
* This is equivalent to
* \f[
\frac{\partial p_C}{\partial \overline{S}_w} =
-\frac{p_e}{\lambda} \overline{S}_w^{-1/\lambda - 1}
\f]
*
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
return twoPhaseSatDPcnw_dSw(params, Sw);
}
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
return - params.entryPressure()/params.lambda() * std::pow(Sw, -1/params.lambda() - 1);
}
/*!
* \brief The partial derivative of the effective saturation with
* regard to the capillary pressure according to Brooks and
* Corey.
*
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar dSw_dpcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
return twoPhaseSatDSw_dpcnw(params, Sw);
}
static Scalar twoPhaseSatDSw_dpcnw(const Params &params, Scalar Sw)
{
assert(pcnw > 0); // required for std::pow
return -params.lambda()/params.entryPressure() * std::pow(pcnw/params.entryPressure(), - params.lambda() - 1);
}
/*!
* \brief The relative permeability for the wetting phase of
* the medium implied by the Brooks-Corey
@ -243,31 +220,25 @@ public:
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
assert(0 <= Sw && Sw <= 1);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
return std::pow(Sw, 2.0/params.lambda() + 3);
const auto& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
/*!
* \brief The derivative of the relative permeability for the
* wetting phase with regard to the wetting saturation of the
* medium implied by the Brooks-Corey parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0 <= Sw && Sw <= 1);
return (2.0/params.lambda() + 3)*std::pow(Sw, 2.0/params.lambda() + 2);
return Toolbox::pow(Sw, 2.0/params.lambda() + 3);
}
/*!
@ -278,39 +249,28 @@ public:
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0 <= Sw && Sw <= 1);
Scalar exponent = 2.0/params.lambda() + 1;
Scalar Sn = 1. - Sw;
return Sn*Sn*(1. - std::pow(Sw, exponent));
const Evaluation Sn = 1. - Sw;
return Sn*Sn*(1. - Toolbox::pow(Sw, exponent));
}
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the medium as implied by the Brooks-Corey
* parameterization.
*
* \param Sw Effective saturation of the wetting phase \f$[-]\f$
* \param params The parameters of the capillary pressure curve
* (for Brooks-Corey: Entry pressure and shape factor)
*/
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{
assert(0 <= Sw && Sw <= 1);
Scalar alpha =
1.0/params.lambda() + 1.0/2 -
Sw*(1.0/params.lambda() + 1.0/2);
return 2.0*(Sw - 1)*(1 + std::pow(Sw, 2.0/params.lambda())*alpha);
}
};
} // namespace Opm

114
opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
View File

@ -26,6 +26,7 @@
#include "EclDefaultMaterialParams.hpp"
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/ErrorMacros.hpp>
@ -133,9 +134,13 @@ public:
const Params &params,
const FluidState &state)
{
values[gasPhaseIdx] = pcgn(params, state);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
values[oilPhaseIdx] = 0;
values[waterPhaseIdx] = - pcnw(params, state);
values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
Valgrind::CheckDefined(values[gasPhaseIdx]);
Valgrind::CheckDefined(values[oilPhaseIdx]);
Valgrind::CheckDefined(values[waterPhaseIdx]);
}
/*!
@ -147,11 +152,13 @@ public:
* p_{c,gn} = p_g - p_n
* \f]
*/
template <class FluidState>
static Scalar pcgn(const Params &params,
const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcgn(const Params &params,
const FluidState &fs)
{
Scalar Sw = 1 - fs.saturation(gasPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = 1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
}
@ -164,12 +171,17 @@ public:
* p_{c,nw} = p_n - p_w
* \f]
*/
template <class FluidState>
static Scalar pcnw(const Params &params,
const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params,
const FluidState &fs)
{
Scalar Sw = fs.saturation(waterPhaseIdx);
return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(waterPhaseIdx));
Valgrind::CheckDefined(Sw);
const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
Valgrind::CheckDefined(result);
return result;
}
/*!
@ -186,9 +198,9 @@ public:
/*!
* \brief The saturation of the gas phase.
*/
template <class FluidState>
static Scalar Sg(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sg(const Params &params,
const FluidState &fluidState)
{
OPM_THROW(std::logic_error, "Not implemented: Sg()");
}
@ -196,9 +208,9 @@ public:
/*!
* \brief The saturation of the non-wetting (i.e., oil) phase.
*/
template <class FluidState>
static Scalar Sn(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params,
const FluidState &fluidState)
{
OPM_THROW(std::logic_error, "Not implemented: Sn()");
}
@ -206,9 +218,9 @@ public:
/*!
* \brief The saturation of the wetting (i.e., water) phase.
*/
template <class FluidState>
static Scalar Sw(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params,
const FluidState &fluidState)
{
OPM_THROW(std::logic_error, "Not implemented: Sw()");
}
@ -233,57 +245,69 @@ public:
const Params &params,
const FluidState &fluidState)
{
values[waterPhaseIdx] = krw(params, fluidState);
values[oilPhaseIdx] = krn(params, fluidState);
values[gasPhaseIdx] = krg(params, fluidState);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
}
/*!
* \brief The relative permeability of the gas phase.
*/
template <class FluidState>
static Scalar krg(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krg(const Params &params,
const FluidState &fluidState)
{
Scalar Sw = 1 - fluidState.saturation(gasPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = 1 - FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
}
/*!
* \brief The relative permeability of the wetting phase.
*/
template <class FluidState>
static Scalar krw(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params,
const FluidState &fluidState)
{
Scalar Sw = fluidState.saturation(waterPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
}
/*!
* \brief The relative permeability of the non-wetting (i.e., oil) phase.
*/
template <class FluidState>
static Scalar krn(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params,
const FluidState &fluidState)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Swco = params.connateWaterSaturation();
Scalar Sw = std::min(1.0, std::max(Swco, fluidState.saturation(waterPhaseIdx)));
//Scalar So = std::min(1.0, std::max(0.0, fluidState.saturation(oilPhaseIdx)));
Scalar Sg = std::min(1.0, std::max(0.0, fluidState.saturation(gasPhaseIdx)));
Evaluation Sw = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(waterPhaseIdx));
// Evaluation So = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(oilPhaseIdx));
Evaluation Sg = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));
Sw = Toolbox::min(1.0, Toolbox::max(Swco, Sw));
//So = Ewoms::Ad::min(1.0, Toolbox::max(0.0, So));
Sg = Toolbox::min(1.0, Toolbox::max(0.0, Sg));
if (Sg + Sw - Swco < 1e-50)
return 1.0; // avoid division by zero
if (Toolbox::value(Sg) + Toolbox::value(Sw) - Swco < 1e-50)
return Toolbox::createConstant(1.0); // avoid division by zero
else {
Scalar kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sg + Sw);
Scalar kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Sw + Swco);
const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sg + Sw);
const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Sw + Swco);
Scalar weightOilWater = (Sw - Swco)/(Sg + Sw - Swco);
Scalar weightGasOil = 1 - weightOilWater;
const auto& weightOilWater = (Sw - Swco)/(Sg + Sw - Swco);
const auto& weightGasOil = 1 - weightOilWater;
Scalar kro = weightOilWater*kro_ow + weightGasOil*kro_go;
return std::min(1.0, std::max(0.0, kro));
Evaluation kro = weightOilWater*kro_ow + weightGasOil*kro_go;
return Toolbox::min(1.0, Toolbox::max(0.0, kro));
}
}
};

150
opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp
View File

@ -122,7 +122,7 @@ public:
phaseIdx));
}
EffLaw::capillaryPressures(values, params, overlayFs);
EffLaw::template capillaryPressures<Container, OverlayFluidState>(values, params, overlayFs);
}
/*!
@ -148,7 +148,7 @@ public:
phaseIdx));
}
EffLaw::relativePermeabilities(values, params, overlayFs);
EffLaw::template relativePermeabilities<Container, OverlayFluidState>(values, params, overlayFs);
}
/*!
@ -162,8 +162,8 @@ public:
* constitutive relation (e.g. Brooks & Corey, van
* Genuchten, linear...)
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<FluidState> OverlayFluidState;
@ -179,14 +179,14 @@ public:
phaseIdx));
}
return EffLaw::pcnw(params, overlayFs);
return EffLaw::template pcnw<OverlayFluidState, Evaluation>(params, overlayFs);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatPcnw(const Params &params, Scalar SwAbs)
template <class Evaluation>
static typename std::enable_if<implementsTwoPhaseSatApi, Evaluation>::type
twoPhaseSatPcnw(const Params &params, const Evaluation& SwAbs)
{
Scalar SwEff = effectiveSaturation(params, SwAbs, Traits::wettingPhaseIdx);
const Evaluation& SwEff = effectiveSaturation(params, SwAbs, Traits::wettingPhaseIdx);
return EffLaw::twoPhaseSatPcnw(params, SwEff);
}
@ -197,7 +197,7 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
EffLaw::saturations(values, params, fs);
EffLaw::template saturations<Container, FluidState>(values, params, fs);
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
values[phaseIdx] = absoluteSaturation(params, values[phaseIdx], phaseIdx);
}
@ -207,27 +207,41 @@ public:
* \brief Calculate wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sw(params, fs), Traits::wettingPhaseIdx); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{
return absoluteSaturation(params,
EffLaw::template Sw<FluidState, Evaluation>(params, fs),
Traits::wettingPhaseIdx);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatSw(const Params &params, Scalar Sw)
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSw(params, Sw), Traits::wettingPhaseIdx); }
template <class Evaluation>
static typename std::enable_if<implementsTwoPhaseSatApi, Evaluation>::type
twoPhaseSatSw(const Params &params, const Evaluation& Sw)
{ return absoluteSaturation(params,
EffLaw::twoPhaseSatSw(params, Sw),
Traits::wettingPhaseIdx); }
/*!
* \brief Calculate non-wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sn(params, fs), Traits::nonWettingPhaseIdx); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{
return absoluteSaturation(params,
EffLaw::template Sn<FluidState, Evaluation>(params, fs),
Traits::nonWettingPhaseIdx);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatSn(const Params &params, Scalar Sw)
{ return absoluteSaturation(params, EffLaw::twoPhaseSatSn(params, Sw), Traits::nonWettingPhaseIdx); }
template <class Evaluation>
static typename std::enable_if<implementsTwoPhaseSatApi, Evaluation>::type
twoPhaseSatSn(const Params &params, const Evaluation& Sw)
{
return absoluteSaturation(params,
EffLaw::twoPhaseSatSn(params, Sw),
Traits::nonWettingPhaseIdx);
}
/*!
* \brief Calculate gas phase saturation given that the rest of
@ -235,51 +249,13 @@ public:
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if< (Traits::numPhases > 2), Evaluation>::type
Sg(const Params &params, const FluidState &fs)
{ return absoluteSaturation(params, EffLaw::Sg(params, fs), Traits::gasPhaseIdx); }
/*!
* \brief Returns the partial derivative of the capillary
* pressure w.r.t the absolute saturation.
*
* In this case the chain rule needs to be applied:
\f[
p_c = p_c( \overline S_w (S_w))
\rightarrow p_c ^\prime = \frac{\partial p_c}{\partial \overline S_w} \frac{\partial \overline S_w}{\partial S_w}
\f]
* \param Sw Absolute saturation of the wetting phase \f$\overline{S}_w\f$.
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of \f$p_c\f$ w.r.t. effective saturation according to EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*/
static Scalar dpC_dSw(const Params &params, Scalar Sw)
{
return EffLaw::dpC_dSw(params, SabsToSeff(params, Sw) )*dSwe_dSw_(params);
}
/*!
* \brief Returns the partial derivative of the absolute
* saturation w.r.t. the capillary pressure.
*
* In this case the chain rule needs to be applied:
\f[
S_w = S_w(\overline{S}_w (p_c) )
\rightarrow S_w^\prime = \frac{\partial S_w}{\partial \overline S_w} \frac{\partial \overline S_w}{\partial p_c}
\f]
*
*
* \param pC Capillary pressure \f$p_C\f$:
* \param params A container object that is populated with the appropriate coefficients for the respective law.
* Therefore, in the (problem specific) spatialParameters first, the material law is chosen, and then the params container
* is constructed accordingly. Afterwards the values are set there, too.
* \return Partial derivative of effective saturation w.r.t. \f$p_c\f$ according to EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*/
static Scalar dSw_dpC(const Params &params, Scalar pC)
{
return EffLaw::dSw_dpC(params, pC)*dSw_dSwe_(params);
return absoluteSaturation(params,
EffLaw::template Sg<FluidState, Evaluation>(params, fs),
Traits::gasPhaseIdx);
}
/*!
@ -291,8 +267,8 @@ public:
* \return Relative permeability of the wetting phase calculated as implied by EffLaw e.g. Brooks & Corey, van Genuchten, linear... .
*
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<FluidState> OverlayFluidState;
@ -308,19 +284,19 @@ public:
phaseIdx));
}
return EffLaw::krw(params, overlayFs);
return EffLaw::template krw<OverlayFluidState, Evaluation>(params, overlayFs);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatKrw(const Params &params, Scalar Sw)
template <class Evaluation>
static typename std::enable_if<implementsTwoPhaseSatApi, Evaluation>::type
twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{ return EffLaw::twoPhaseSatKrw(params, effectiveSaturation(params, Sw, Traits::nonWettingPhaseIdx)); }
/*!
* \brief The relative permeability of the non-wetting phase.
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<FluidState> OverlayFluidState;
@ -336,12 +312,12 @@ public:
phaseIdx));
}
return EffLaw::krn(params, overlayFs);
return EffLaw::template krn<OverlayFluidState, Evaluation>(params, overlayFs);
}
template <class ScalarT = Scalar>
static typename std::enable_if<implementsTwoPhaseSatApi, ScalarT>::type
twoPhaseSatKrn(const Params &params, Scalar Sw)
template <class Evaluation>
static typename std::enable_if<implementsTwoPhaseSatApi, Evaluation>::type
twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{ return EffLaw::twoPhaseSatKrn(params, effectiveSaturation(params, Sw, Traits::nonWettingPhaseIdx)); }
/*!
@ -349,8 +325,8 @@ public:
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if< (Traits::numPhases > 2), Evaluation>::type
krg(const Params &params, const FluidState &fs)
{
typedef Opm::SaturationOverlayFluidState<FluidState> OverlayFluidState;
@ -367,20 +343,22 @@ public:
phaseIdx));
}
return EffLaw::krg(params, overlayFs);
return EffLaw::template krg<OverlayFluidState, Evaluation>(params, overlayFs);
}
/*!
* \brief Convert an absolute saturation to an effective one.
*/
static Scalar effectiveSaturation(const Params &params, Scalar S, int phaseIdx)
{ return (S - params.residualSaturation(phaseIdx))/(1 - params.sumResidualSaturations()); }
template <class Evaluation>
static Evaluation effectiveSaturation(const Params &params, const Evaluation& S, int phaseIdx)
{ return (S - params.residualSaturation(phaseIdx))/(1.0 - params.sumResidualSaturations()); }
/*!
* \brief Convert an effective saturation to an absolute one.
*/
static Scalar absoluteSaturation(const Params &params, Scalar S, int phaseIdx)
{ return S*(1 - params.sumResidualSaturations()) + params.residualSaturation(phaseIdx); }
template <class Evaluation>
static Evaluation absoluteSaturation(const Params &params, const Evaluation& S, int phaseIdx)
{ return S*(1.0 - params.sumResidualSaturations()) + params.residualSaturation(phaseIdx); }
private:
/*!

169
opm/material/fluidmatrixinteractions/LinearMaterial.hpp
View File

@ -25,8 +25,8 @@
#include "LinearMaterialParams.hpp"
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/ErrorMacros.hpp>
@ -97,8 +97,12 @@ public:
const Params &params,
const FluidState &state)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
for (int phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
Scalar S = state.saturation(phaseIdx);
const Evaluation& S =
FsToolbox::template toLhs<Evaluation>(state.saturation(phaseIdx));
Valgrind::CheckDefined(S);
values[phaseIdx] =
@ -126,47 +130,55 @@ public:
const Params &params,
const FluidState &state)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
for (int phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) {
Scalar S = state.saturation(phaseIdx);
const Evaluation& S =
FsToolbox::template toLhs<Evaluation>(state.saturation(phaseIdx));
Valgrind::CheckDefined(S);
values[phaseIdx] = std::max(std::min(S,1.0),0.0);
values[phaseIdx] = Toolbox::max(Toolbox::min(S,1.0),0.0);
}
}
/*!
* \brief The difference between the pressures of the non-wetting and wetting phase.
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
Scalar S = fs.saturation(Traits::wettingPhaseIdx);
Valgrind::CheckDefined(S);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
Valgrind::CheckDefined(Sw);
Scalar wPhasePressure =
S*params.pcMaxSat(Traits::wettingPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::wettingPhaseIdx);
const Evaluation& wPhasePressure =
Sw*params.pcMaxSat(Traits::wettingPhaseIdx) +
(1.0 - Sw)*params.pcMinSat(Traits::wettingPhaseIdx);
S = fs.saturation(Traits::nonWettingPhaseIdx);
Valgrind::CheckDefined(S);
const Evaluation& Sn =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
Valgrind::CheckDefined(Sn);
Scalar nPhasePressure =
S*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nonWettingPhaseIdx);
const Evaluation& nPhasePressure =
Sn*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
(1.0 - Sn)*params.pcMinSat(Traits::nonWettingPhaseIdx);
return nPhasePressure - wPhasePressure;
}
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatPcnw(const Params &params, Scalar Sw)
template <class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{
Scalar wPhasePressure =
const Evaluation& wPhasePressure =
Sw*params.pcMaxSat(Traits::wettingPhaseIdx) +
(1.0 - Sw)*params.pcMinSat(Traits::wettingPhaseIdx);
Scalar nPhasePressure =
const Evaluation& nPhasePressure =
(1.0 - Sw)*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
Sw*params.pcMinSat(Traits::nonWettingPhaseIdx);
@ -177,26 +189,26 @@ public:
* \brief Calculate wetting phase saturation given that the rest
* of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{ OPM_THROW(std::runtime_error, "Not implemented: Sw()"); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatSw(const Params &params, Scalar Sw)
template <class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
twoPhaseSatSw(const Params &params, const Evaluation& Sw)
{ OPM_THROW(std::runtime_error, "Not implemented: twoPhaseSatSw()"); }
/*!
* \brief Calculate non-wetting liquid phase saturation given that
* the rest of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{ OPM_THROW(std::runtime_error, "Not implemented: Sn()"); }
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatSn(const Params &params, Scalar Sw)
template <class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
twoPhaseSatSn(const Params &params, const Evaluation& Sw)
{ OPM_THROW(std::runtime_error, "Not implemented: twoPhaseSatSn()"); }
/*!
@ -205,67 +217,100 @@ public:
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
Sg(const Params &params, const FluidState &fs)
{ OPM_THROW(std::runtime_error, "Not implemented: Sg()"); }
/*!
* \brief The relative permability of the wetting phase
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::wettingPhaseIdx))); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, Sw)); }
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
}
template <class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
}
/*!
* \brief The relative permability of the liquid non-wetting phase
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::nonWettingPhaseIdx))); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
template <class ScalarT = Scalar>
static typename std::enable_if<Traits::numPhases == 2, ScalarT>::type
twoPhaseSatKrn(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, 1 - Sw)); }
const Evaluation& Sn =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sn));
}
template <class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 2, Evaluation>::type
twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
}
/*!
* \brief The relative permability of the gas phase
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
krg(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::gasPhaseIdx))); }
{
typedef MathToolbox<Evaluation> Toolbox;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sg =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sg));
}
/*!
* \brief The difference between the pressures of the gas and the non-wetting phase.
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
template <class FluidState, class Evaluation = Scalar>
static typename std::enable_if<Traits::numPhases == 3, Evaluation>::type
pcgn(const Params &params, const FluidState &fs)
{
Scalar S = fs.saturation(Traits::nonWettingPhaseIdx);
Valgrind::CheckDefined(S);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar nPhasePressure =
S*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::nonWettingPhaseIdx);
const Evaluation& Sn =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
Valgrind::CheckDefined(Sn);
S = fs.saturation(Traits::gasPhaseIdx);
Valgrind::CheckDefined(S);
const Evaluation& nPhasePressure =
Sn*params.pcMaxSat(Traits::nonWettingPhaseIdx) +
(1.0 - Sn)*params.pcMinSat(Traits::nonWettingPhaseIdx);
Scalar gPhasePressure =
S*params.pcMaxSat(Traits::gasPhaseIdx) +
(1.0 - S)*params.pcMinSat(Traits::gasPhaseIdx);
const Evaluation& Sg =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::gasPhaseIdx));
Valgrind::CheckDefined(Sg);
const Evaluation& gPhasePressure =
Sg*params.pcMaxSat(Traits::gasPhaseIdx) +
(1.0 - Sg)*params.pcMinSat(Traits::gasPhaseIdx);
return gPhasePressure - nPhasePressure;
}

19
opm/material/fluidmatrixinteractions/MaterialTraits.hpp
View File

@ -35,13 +35,16 @@ namespace Opm {
*
* This traits class is intended to be used by the NullMaterial
*/
template <class ScalarT, int numPhasesV>
template <class ScalarT, int numPhasesV, class EvaluationT = ScalarT>
class NullMaterialTraits
{
public:
//! The type used for scalar floating point values
typedef ScalarT Scalar;
//! Representation of a function evaluation and all its relevant derivatives
typedef EvaluationT Evaluation;
//! The number of fluid phases
static const int numPhases = numPhasesV;
};
@ -51,21 +54,24 @@ public:
*
* \brief A generic traits class for two-phase material laws.
*/
template <class ScalarT, int wettinPhaseIdxV, int nonWettingasPhaseIdxV>
template <class ScalarT, int wettingPhaseIdxV, int nonWettingPhaseIdxV, class EvaluationT = ScalarT>
class TwoPhaseMaterialTraits
{
public:
//! The type used for scalar floating point values
typedef ScalarT Scalar;
//! Representation of a function evaluation and all its relevant derivatives
typedef EvaluationT Evaluation;
//! The number of fluid phases
static const int numPhases = 2;
//! The index of the wetting phase
static const int wettingPhaseIdx = wettinPhaseIdxV;
static const int wettingPhaseIdx = wettingPhaseIdxV;
//! The index of the non-wetting phase
static const int nonWettingPhaseIdx = nonWettingasPhaseIdxV;
static const int nonWettingPhaseIdx = nonWettingPhaseIdxV;
// some safety checks...
static_assert(wettingPhaseIdx != nonWettingPhaseIdx,
@ -77,13 +83,16 @@ public:
*
* \brief A generic traits class for three-phase material laws.
*/
template <class ScalarT, int wettingPhaseIdxV, int nonWettingasPhaseIdxV, int gasPhaseIdxV>
template <class ScalarT, int wettingPhaseIdxV, int nonWettingasPhaseIdxV, int gasPhaseIdxV, class EvaluationT = ScalarT>
class ThreePhaseMaterialTraits
{
public:
//! The type used for scalar floating point values
typedef ScalarT Scalar;
//! Representation of a function evaluation and all its relevant derivatives
typedef EvaluationT Evaluation;
//! The number of fluid phases
static const int numPhases = 3;

127
opm/material/fluidmatrixinteractions/NullMaterial.hpp
View File

@ -27,6 +27,7 @@
#include <opm/material/common/ErrorMacros.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <algorithm>
@ -107,50 +108,56 @@ public:
template <class ContainerT, class FluidState>
static void relativePermeabilities(ContainerT &values,
const Params &params,
const FluidState &state)
const FluidState &fluidState)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
values[phaseIdx] = std::max(std::min(state.saturation(phaseIdx),1.0),0.0);
const Evaluation& S =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(phaseIdx));
values[phaseIdx] = Toolbox::max(Toolbox::min(S, 1.0), 0.0);
}
}
/*!
* \brief The difference between the pressures of the non-wetting and wetting phase.
*/
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if<(numPhases > 1), ScalarT>::type
pcnw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if<(numPhases > 1), Evaluation>::type
pcnw(const Params &params, const FluidState &fluidState)
{ return 0; }
template <class ScalarT = Scalar>
static typename std::enable_if<numPhases == 2, ScalarT>::type
twoPhaseSatPcnw(const Params &params, Scalar Sw)
template <class Evaluation>
static typename std::enable_if<numPhases == 2, Evaluation>::type
twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{ return 0; }
/*!
* \brief Calculate wetting phase saturation given that the rest
* of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Scalar Sw(const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not defined: Sw()"); }
template <class ScalarT = Scalar>
static typename std::enable_if<numPhases == 2, ScalarT>::type
twoPhaseSatSw(const Params &params, Scalar pcnw)
template <class Evaluation>
static typename std::enable_if<numPhases == 2, Evaluation>::type
twoPhaseSatSw(const Params &params, const Evaluation& pcnw)
{ OPM_THROW(std::logic_error, "Not defined: twoPhaseSatSw()"); }
/*!
* \brief Calculate non-wetting phase saturation given that the
* rest of the fluid state has been initialized
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Scalar Sn(const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not defined: Sn()"); }
template <class ScalarT = Scalar>
static typename std::enable_if<numPhases == 2, ScalarT>::type
twoPhaseSatSn(const Params &params, Scalar pcnw)
template <class Evaluation>
static typename std::enable_if<numPhases == 2, Evaluation>::type
twoPhaseSatSn(const Params &params, const Evaluation& pcnw)
{ OPM_THROW(std::logic_error, "Not defined: twoPhaseSatSn()"); }
/*!
@ -159,55 +166,87 @@ public:
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if< (numPhases > 2), ScalarT>::type
Sg(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if< (numPhases > 2), Evaluation>::type
Sg(const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not defined: Sg()"); }
/*!
* \brief The relative permability of the wetting phase
*/
template <class FluidState, class ScalarT = Scalar>
static typename std::enable_if<(numPhases > 1), ScalarT>::type
krw(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::wettingPhaseIdx))); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if<(numPhases > 1), Evaluation>::type
krw(const Params &params, const FluidState &fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
template <class ScalarT = Scalar>
static typename std::enable_if<numPhases == 2, ScalarT>::type
twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, Sw)); }
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
}
template <class Evaluation>
static typename std::enable_if<numPhases == 2, Evaluation>::type
twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, Sw));
}
/*!
* \brief The relative permability of the liquid non-wetting phase
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if<(numPhases > 1), ScalarT>::type
krn(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::nonWettingPhaseIdx))); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if<(numPhases > 1), Evaluation>::type
krn(const Params &params, const FluidState &fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
template <class ScalarT = Scalar>
static typename std::enable_if<numPhases == 2, ScalarT>::type
twoPhaseSatKrn(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, 1 - Sw)); }
const Evaluation& Sn =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(Traits::nonWettingPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sn));
}
template <class Evaluation>
static typename std::enable_if<numPhases == 2, Evaluation>::type
twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, 1.0 - Sw));
}
/*!
* \brief The relative permability of the gas phase
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (numPhases > 2), ScalarT>::type
krg(const Params &params, const FluidState &fs)
{ return std::max(0.0, std::min(1.0, fs.saturation(Traits::gasPhaseIdx))); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if< (numPhases > 2), Evaluation>::type
krg(const Params &params, const FluidState &fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sg =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(Traits::gasPhaseIdx));
return Toolbox::max(0.0, Toolbox::min(1.0, Sg));
}
/*!
* \brief The difference between the pressures of the gas and the non-wetting phase.
*
* This method is only available for at least three fluid phases
*/
template <class FluidState, class ScalarT=Scalar>
static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type
pcgn(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static typename std::enable_if< (Traits::numPhases > 2), Evaluation>::type
pcgn(const Params &params, const FluidState &fluidState)
{ return 0; }
};
} // namespace Opm

132
opm/material/fluidmatrixinteractions/ParkerLenhard.hpp
View File

@ -299,7 +299,9 @@ public:
template <class FluidState>
static void update(Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar Sw = FsToolbox::value(fs.saturation(Traits::wettingPhaseIdx));
if (Sw > 1 - 1e-5) {
// if the absolute saturation is almost 1,
@ -315,7 +317,7 @@ public:
// calculate the apparent saturation on the MIC and MDC
// which yield the same capillary pressure as the
// Sw at the current scanning curve
Scalar pc = pcnw(params, fs);
Scalar pc = pcnw<FluidState, Scalar>(params, fs);
Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc);
Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc);
@ -339,8 +341,10 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
}
/*!
@ -358,25 +362,37 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
}
/*!
* \brief Returns the capillary pressure dependend on
* the phase saturations.
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
// calculate the current apparent saturation
ScanningCurve *sc = findScanningCurve_(params, Sw);
ScanningCurve *sc = findScanningCurve_(params, Toolbox::value(Sw));
// calculate the apparant saturation
Scalar Sw_app = absoluteToApparentSw_(params, Sw);
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
// if the apparent saturation exceeds the 'legal' limits,
// we also the underlying material law decide what to do.
@ -388,88 +404,98 @@ public:
// drainage curve
Scalar SwAppCurSC = absoluteToApparentSw_(params, sc->Sw());
Scalar SwAppPrevSC = absoluteToApparentSw_(params, sc->prev()->Sw());
Scalar pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
const Evaluation& pos = (Sw_app - SwAppCurSC)/(SwAppPrevSC - SwAppCurSC);
if (sc->isImbib()) {
Scalar SwMic =
const Evaluation& SwMic =
pos * (sc->prev()->SwMic() - sc->SwMic()) + sc->SwMic();
return VanGenuchten::twoPhaseSatPcnw(params.micParams(), SwMic);
}
else { // sc->isDrain()
Scalar SwMdc =
const Evaluation& SwMdc =
pos*(sc->prev()->SwMdc() - sc->SwMdc()) + sc->SwMdc();
return VanGenuchten::twoPhaseSatPcnw(params.mdcParams(), SwMdc);
}
}
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{ OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatDPcnw_dSw()"); }
/*!
* \brief Calculate the wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{ OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::Sw()"); }
static Scalar twoPhaseSatSw(const Params &params, Scalar pc)
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pc)
{ OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::twoPhaseSatSw()"); }
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw<FluidState, Evaluation>(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pc)
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pc)
{ OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::twoPhaseSatSn()"); }
/*!
* \brief The relative permeability for the wetting phase of
* the medium.
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
// for the effective permeability we only use Land's law and
// the relative permeability of the main drainage curve.
Scalar Sw_app = absoluteToApparentSw_(params, Sw);
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
return VanGenuchten::twoPhaseSatKrw(params.mdcParams(), Sw_app);
}
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
{ OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatDKrw_dSw()"); }
/*!
* \brief The relative permeability for the non-wetting phase
* of the params.
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, fs.saturation(Traits::wettingPhaseIdx)); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
// for the effective permeability we only use Land's law and
// the relative permeability of the main drainage curve.
Scalar Sw_app = absoluteToApparentSw_(params, Sw);
const Evaluation& Sw_app = absoluteToApparentSw_(params, Sw);
return VanGenuchten::twoPhaseSatKrn(params.mdcParams(), Sw_app);
}
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{ OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatDKrn_dSw()"); }
/*!
* \brief Convert an absolute wetting saturation to an apparent one.
*/
static Scalar absoluteToApparentSw_(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation absoluteToApparentSw_(const Params &params, const Evaluation& Sw)
{
return effectiveToApparentSw_(params, absoluteToEffectiveSw_(params, Sw));
}
@ -484,7 +510,8 @@ private:
* is constructed accordingly. Afterwards the values are set there, too.
* \return Effective saturation of the wetting phase.
*/
static Scalar absoluteToEffectiveSw_(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation absoluteToEffectiveSw_(const Params &params, const Evaluation& Sw)
{ return (Sw - params.SwrPc())/(1 - params.SwrPc()); }
/*!
@ -496,12 +523,14 @@ private:
* is constructed accordingly. Afterwards the values are set there, too.
* \return Absolute saturation of the non-wetting phase.
*/
static Scalar effectiveToAbsoluteSw_(const Params &params, Scalar Swe)
template <class Evaluation>
static Evaluation effectiveToAbsoluteSw_(const Params &params, const Evaluation& Swe)
{ return Swe*(1 - params.SwrPc()) + params.SwrPc(); }
// return the effctive residual non-wetting saturation, given an
// effective wetting saturation
static Scalar computeCurrentSnr_(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation computeCurrentSnr_(const Params &params, const Evaluation& Sw)
{
// regularize
if (Sw > 1 - params.Snr())
@ -514,7 +543,7 @@ private:
// use Land's law
Scalar R = 1.0/params.Snr() - 1;
Scalar curSnr = (1 - Sw)/(1 + R*(1 - Sw));
const Evaluation& curSnr = (1 - Sw)/(1 + R*(1 - Sw));
// the current effective residual non-wetting saturation must
// be smaller than the residual non-wetting saturation
@ -525,9 +554,10 @@ private:
// returns the trapped effective non-wetting saturation for a
// given wetting phase saturation
static Scalar trappedEffectiveSn_(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation trappedEffectiveSn_(const Params &params, const Evaluation& Sw)
{
Scalar Swe = absoluteToEffectiveSw_(params, Sw);
const Evaluation& Swe = absoluteToEffectiveSw_(params, Sw);
Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
Scalar Snre = absoluteToEffectiveSw_(params, params.currentSnr());
@ -536,7 +566,8 @@ private:
// returns the apparent saturation of the wetting phase depending
// on the effective saturation
static Scalar effectiveToApparentSw_(const Params &params, Scalar Swe)
template <class Evaluation>
static Evaluation effectiveToApparentSw_(const Params &params, const Evaluation& Swe)
{
if (params.pisc() == NULL ||
Swe <= absoluteToEffectiveSw_(params, params.pisc()->Sw()))
@ -554,7 +585,8 @@ private:
}
// Returns the effective saturation to a given apparent one
static Scalar apparentToEffectiveSw_(const Params &params, Scalar Swapp)
template <class Evaluation>
static Evaluation apparentToEffectiveSw_(const Params &params, const Evaluation& Swapp)
{
Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw());
if (params.pisc() == NULL || Swapp <= SwePisc) {

168
opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp
View File

@ -27,6 +27,7 @@
#include <opm/material/common/ErrorMacros.hpp>
#include <opm/material/common/Exceptions.hpp>
#include <opm/material/common/MathToolbox.hpp>
#include <algorithm>
#include <cmath>
@ -94,8 +95,10 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
}
/*!
@ -112,17 +115,21 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
}
/*!
* \brief The capillary pressure-saturation curve
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
@ -130,117 +137,94 @@ public:
/*!
* \brief The saturation-capillary pressure curve
*/
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{ return eval_(params.SwSamples(), params.pcnwSamples(), Sw); }
/*!
* \brief The saturation-capillary pressure curve
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{ OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
static Scalar twoPhaseSatSw(const Params &params, Scalar pC)
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pC)
{ OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); }
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw<FluidState, Scalar>(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pC)
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pC)
{ return 1 - twoPhaseSatSw(params, pC); }
/*!
* \brief The partial derivative of the capillary pressure with
* regard to the saturation
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
assert(0 < Sw && Sw < 1);
return evalDeriv_(params.SwSamples(),
params.pcnwSamples(),
Sw);
}
/*!
* \brief The relative permeability for the wetting phase of the
* porous medium
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, eval_(params.SwSamples(), params.krwSamples(), Sw))); }
/*!
* \brief The derivative of the relative permeability of the
* wetting phase in regard to the wetting saturation of the
* porous medium
*/
template <class FluidState>
static Scalar dKrw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
return evalDeriv_(params.SwSamples(),
params.krwSamples(),
Sw);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
const auto& res = eval_(params.SwSamples(), params.krwSamples(), Sw);
return Toolbox::max(0.0, Toolbox::min(1.0, res));
}
/*!
* \brief The relative permeability for the non-wetting phase
* of the porous medium
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
return std::max(0.0, std::min(1.0, eval_(params.SwSamples(),
params.krnSamples(),
Sw)));
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the porous medium
*/
template <class FluidState>
static Scalar dKrn_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
if (Sw < params.SwSamples().front() || Sw > params.SwSamples().back())
return 0.0;
return evalDeriv_(params.SwSamples(),
params.krnSamples(),
Sw);
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, eval_(params.SwSamples(),
params.krnSamples(),
Sw)));
}
private:
static Scalar eval_(const ValueVector &xValues,
const ValueVector &yValues,
Scalar x)
template <class Evaluation>
static Evaluation eval_(const ValueVector &xValues,
const ValueVector &yValues,
const Evaluation& x)
{
if (x < xValues.front())
return yValues.front();
if (x > xValues.back())
return yValues.back();
typedef MathToolbox<Evaluation> Toolbox;
int segIdx = findSegmentIndex_(xValues, x);
if (Toolbox::value(x) < xValues.front())
return Toolbox::createConstant(yValues.front());
if (Toolbox::value(x) > xValues.back())
return Toolbox::createConstant(yValues.back());
int segIdx = findSegmentIndex_(xValues, Toolbox::value(x));
Scalar x0 = xValues[segIdx];
Scalar x1 = xValues[segIdx + 1];
@ -248,14 +232,24 @@ private:
Scalar y0 = yValues[segIdx];
Scalar y1 = yValues[segIdx + 1];
return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
Scalar m = (y1 - y0)/(x1 - x0);
return y0 + (x - x0)*m;
}
static Scalar evalDeriv_(const ValueVector &xValues,
const ValueVector &yValues,
Scalar x)
template <class Evaluation>
static Evaluation evalDeriv_(const ValueVector &xValues,
const ValueVector &yValues,
const Evaluation& x)
{
int segIdx = findSegmentIndex_(xValues, x);
typedef MathToolbox<Evaluation> Toolbox;
if (Toolbox::value(x) < xValues.front())
return Toolbox::createConstant(0.0);
if (Toolbox::value(x) > xValues.back())
return Toolbox::createConstant(0.0);
int segIdx = findSegmentIndex_(xValues, Toolbox::value(x));
Scalar x0 = xValues[segIdx];
Scalar x1 = xValues[segIdx + 1];
@ -263,7 +257,7 @@ private:
Scalar y0 = yValues[segIdx];
Scalar y1 = yValues[segIdx + 1];
return (y1 - y0)/(x1 - x0);
return Toolbox::createConstant((y1 - y0)/(x1 - x0));
}
static int findSegmentIndex_(const ValueVector &xValues, Scalar x)

42
opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp
View File

@ -97,8 +97,18 @@ public:
*
* This curve is assumed to depend on the wetting phase saturation
*/
void setPcnwSamples(const ValueVector& SwValues, const ValueVector& values)
{ SwSamples_ = SwValues; pcwnSamples_ = values; }
template <class Container>
void setPcnwSamples(const Container& SwValues, const Container& values)
{
assert(SwValues.size() == values.size());
int n = SwValues.size();
SwSamples_.resize(n);
pcwnSamples_.resize(n);
std::copy(SwValues.begin(), SwValues.end(), SwSamples_.begin());
std::copy(values.begin(), values.end(), pcwnSamples_.begin());
}
/*!
* \brief Return the sampling points for the relative permeability
@ -115,8 +125,18 @@ public:
*
* This curve is assumed to depend on the wetting phase saturation
*/
void setKrwSamples(const ValueVector& SwValues, const ValueVector& values)
{ SwSamples_ = SwValues; krwSamples_ = values; }
template <class Container>
void setKrwSamples(const Container& SwValues, const Container& values)
{
assert(SwValues.size() == values.size());
int n = SwValues.size();
SwSamples_.resize(n);
krwSamples_.resize(n);
std::copy(SwValues.begin(), SwValues.end(), SwSamples_.begin());
std::copy(values.begin(), values.end(), krwSamples_.begin());
}
/*!
* \brief Return the sampling points for the relative permeability
@ -133,8 +153,18 @@ public:
*
* This curve is assumed to depend on the wetting phase saturation
*/
void setKrnSamples(const ValueVector& SwValues, const ValueVector& values)
{ SwSamples_ = SwValues; krnSamples_ = values; }
template <class Container>
void setKrnSamples(const Container& SwValues, const Container& values)
{
assert(SwValues.size() == values.size());
int n = SwValues.size();
SwSamples_.resize(n);
krnSamples_.resize(n);
std::copy(SwValues.begin(), SwValues.end(), SwSamples_.begin());
std::copy(values.begin(), values.end(), krnSamples_.begin());
}
private:
#ifndef NDEBUG

192
opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp
View File

@ -115,8 +115,10 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
}
/*!
@ -126,7 +128,9 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = Sw(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
@ -143,8 +147,10 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
}
/*!
@ -171,13 +177,19 @@ public:
*
* \sa BrooksCorey::pcnw
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
const Scalar Sthres = params.thresholdSw();
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{
const Scalar Sthres = params.pcnwLowSw();
// make sure that the capilary pressure observes a
// derivative != 0 for 'illegal' saturations. This is
@ -186,13 +198,13 @@ public:
// saturation moving to the right direction if it
// temporarily is in an 'illegal' range.
if (Sw <= Sthres) {
Scalar m = BrooksCorey::twoPhaseSatDPcnw_dSw(params, Sthres);
Scalar pcnw_SwLow = BrooksCorey::twoPhaseSatPcnw(params, Sthres);
Scalar m = params.pcnwSlopeLow();
Scalar pcnw_SwLow = params.pcnwLow();
return pcnw_SwLow + m*(Sw - Sthres);
}
else if (Sw > 1.0) {
Scalar m = BrooksCorey::twoPhaseSatDPcnw_dSw(params, 1.0);
Scalar pcnw_SwHigh = params.entryPressure();
else if (Sw >= 1.0) {
Scalar m = params.pcnwSlopeHigh();
Scalar pcnw_SwHigh = params.pcnwHigh();
return pcnw_SwHigh + m*(Sw - 1.0);
}
@ -207,23 +219,28 @@ public:
*
* This is the inverse of the pcnw() method.
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{
Scalar pcnw = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
return twoPhaseSatSw(params, pcnw);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& pC =
FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pcnw)
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pcnw)
{
const Scalar Sthres = params.thresholdSw();
const Scalar Sthres = params.pcnwLowSw();
// calculate the saturation which corrosponds to the
// saturation in the non-regularized version of the
// Brooks-Corey law. If the input capillary pressure is
// smaller than the entry pressure, make sure that we will
// regularize.
Scalar Sw = 1.5;
Evaluation Sw = 1.5;
if (pcnw >= params.entryPressure())
Sw = BrooksCorey::twoPhaseSatSw(params, pcnw);
@ -235,13 +252,13 @@ public:
// temporarily is in an 'illegal' range.
if (Sw <= Sthres) {
// invert the low saturation regularization of pcnw()
Scalar m = BrooksCorey::twoPhaseSatDPcnw_dSw(params, Sthres);
Scalar pcnw_SwLow = BrooksCorey::twoPhaseSatPcnw(params, Sthres);
Scalar m = params.pcnwSlopeLow();
Scalar pcnw_SwLow = params.pcnwLow();
return Sthres + (pcnw - pcnw_SwLow)/m;
}
else if (Sw > 1.0) {
Scalar m = BrooksCorey::twoPhaseSatDPcnw_dSw(params, 1.0);
Scalar pcnw_SwHigh = BrooksCorey::twoPhaseSatPcnw(params, 1.0);
Scalar m = params.pcnwSlopeHigh();
Scalar pcnw_SwHigh = params.pcnwHigh();
return 1.0 + (pcnw - pcnw_SwHigh)/m;;
}
@ -252,67 +269,14 @@ public:
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw<FluidState, Evaluation>(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pcnw)
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pcnw)
{ return 1 - twoPhaseSatSw(params, pcnw); }
/*!
* \brief The derivative of the regularized Brooks-Corey capillary
* pressure-saturation curve.
*/
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
const Scalar Sthres = params.thresholdSw();
// derivative of the regualarization
if (Sw <= Sthres) {
// calculate the slope of the straight line used in pcnw()
Scalar m = BrooksCorey::twoPhaseSatDPcnw_dSw(params, Sthres);
return m;
}
else if (Sw > 1.0) {
// calculate the slope of the straight line used in pcnw()
Scalar m = BrooksCorey::twoPhaseSatDPcnw_dSw(params, 1.0);
return m;
}
return BrooksCorey::twoPhaseSatDPcnw_dSw(params, Sw);
}
/*!
* \brief The derivative of the regularized Brooks-Corey
* saturation-capillary pressure curve.
*/
static Scalar twoPhaseSatDSw_dpcnw(const Params &params, Scalar pcnw)
{
const Scalar Sthres = params.thresholdSw();
// calculate the saturation which corresponds to the
// saturation in the non-regularized version of the
// Brooks-Corey law
Scalar Sw;
if (pcnw < params.entryPressure())
Sw = 1.5; // make sure we regularize (see below)
else
Sw = BrooksCorey::Sw(params, pcnw);
// derivative of the regularization
if (Sw <= Sthres) {
// calculate the slope of the straight line used in pcnw()
Scalar m = BrooksCorey::dPcnw_dSw(params, Sthres);
return 1/m;
}
else if (Sw > 1.0) {
// calculate the slope of the straight line used in pcnw()
Scalar m = BrooksCorey::dPcnw_dSw(params, 1.0);
return 1/m;
}
return 1.0/BrooksCorey::dPcnw_dSw(params, Sw);
}
/*!
* \brief Regularized version of the relative permeability of the
* wetting phase of the Brooks-Corey curves.
@ -327,26 +291,26 @@ public:
*
* \sa BrooksCorey::krw
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
if (Sw <= 0.0)
return 0.0;
else if (Sw >= 1.0)
return 1.0;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
return BrooksCorey::twoPhaseSatKrw(params, Sw);
const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
if (Sw <= 0.0 || Sw >= 1.0)
return 0.0;
typedef MathToolbox<Evaluation> Toolbox;
return BrooksCorey::twoPhaseSatDKrw_dSw(params, Sw);
if (Sw <= 0.0)
return Toolbox::createConstant(0.0);
else if (Sw >= 1.0)
return Toolbox::createConstant(1.0);
return BrooksCorey::twoPhaseSatKrw(params, Sw);
}
/*!
@ -363,28 +327,28 @@ public:
*
* \sa BrooksCorey::krn
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
if (Sw >= 1.0)
return 0.0;
return Toolbox::createConstant(0.0);
else if (Sw <= 0.0)
return 1.0;
return Toolbox::createConstant(1.0);
return BrooksCorey::twoPhaseSatKrn(params, Sw);
}
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{
if (Sw <= 0.0 || Sw >= 1.0)
return 0.0;
return BrooksCorey::twoPhaseSatDKrn_dSw(params, Sw);
}
};
} // namespace Opm

96
opm/material/fluidmatrixinteractions/RegularizedBrooksCoreyParams.hpp
View File

@ -23,8 +23,11 @@
#ifndef OPM_REGULARIZED_BROOKS_COREY_PARAMS_HPP
#define OPM_REGULARIZED_BROOKS_COREY_PARAMS_HPP
#include "BrooksCorey.hpp"
#include "BrooksCoreyParams.hpp"
#include <dune/common/deprecated.hh>
#include <cassert>
namespace Opm {
@ -38,6 +41,7 @@ template <class TraitsT>
class RegularizedBrooksCoreyParams : public Opm::BrooksCoreyParams<TraitsT>
{
typedef Opm::BrooksCoreyParams<TraitsT> BrooksCoreyParams;
typedef Opm::BrooksCorey<TraitsT, RegularizedBrooksCoreyParams> BrooksCorey;
typedef typename TraitsT::Scalar Scalar;
public:
@ -45,7 +49,7 @@ public:
RegularizedBrooksCoreyParams()
: BrooksCoreyParams()
, SwThres_(1e-2)
, pcnwLowSw_(0.01)
{
#ifndef NDEBUG
finalized_ = false;
@ -54,7 +58,7 @@ public:
RegularizedBrooksCoreyParams(Scalar entryPressure, Scalar lambda)
: BrooksCoreyParams(entryPressure, lambda)
, SwThres_(1e-2)
, pcnwLowSw_(0.01)
{ finalize(); }
/*!
@ -64,24 +68,65 @@ public:
void finalize()
{
BrooksCoreyParams::finalize();
pcnwLow_ = BrooksCorey::twoPhaseSatPcnw(*this, pcnwLowSw_);
pcnwSlopeLow_ = dPcnw_dSw_(pcnwLowSw_);
pcnwHigh_ = BrooksCorey::twoPhaseSatPcnw(*this, 1.0);
pcnwSlopeHigh_ = dPcnw_dSw_(1.0);
#ifndef NDEBUG
finalized_ = true;
#endif
}
/*!
* \brief Return the threshold saturation below which the capillary pressure
* is regularized.
* \brief Return the threshold saturation below which the
* capillary pressure is regularized.
*/
Scalar thresholdSw() const
{ assertFinalized_(); return SwThres_; }
Scalar pcnwLowSw() const
{ assertFinalized_(); return pcnwLowSw_; }
/*!
* \brief Set the threshold saturation below which the capillary pressure
* is regularized.
* \brief Return the capillary pressure at the low threshold
* saturation of the wetting phase.
*/
void setThresholdSw(Scalar value)
{ SwThres_ = value; }
Scalar pcnwLow() const
{ assertFinalized_(); return pcnwLow_; }
/*!
* \brief Return the slope capillary pressure curve if Sw is
* smaller or equal to the low threshold saturation.
*
* For this case, we extrapolate the curve using a straight line.
*/
Scalar pcnwSlopeLow() const
{ assertFinalized_(); return pcnwSlopeLow_; }
/*!
* \brief Set the threshold saturation below which the capillary
* pressure is regularized.
*/
void setPcLowSw(Scalar value)
{ pcnwLowSw_ = value; }
void setThresholdSw(Scalar value) DUNE_DEPRECATED_MSG("this method has been renamed to setPcLowSw()")
{ pcnwLowSw_ = value; }
/*!
* \brief Return the capillary pressure at the high threshold
* saturation of the wetting phase.
*/
Scalar pcnwHigh() const
{ assertFinalized_(); return pcnwHigh_; }
/*!
* \brief Return the slope capillary pressure curve if Sw is
* larger or equal to 1.
*
* For this case, we extrapolate the curve using a straight line.
*/
Scalar pcnwSlopeHigh() const
{ assertFinalized_(); return pcnwSlopeHigh_; }
private:
#ifndef NDEBUG
@ -94,7 +139,36 @@ private:
{ }
#endif
Scalar SwThres_;
Scalar dPcnw_dSw_(Scalar Sw) const
{
// use finite differences to calculate the derivative w.r.t. Sw of the
// unregularized curve's capillary pressure.
const Scalar eps = 1e-7;
Scalar delta = 0.0;
Scalar pc1;
Scalar pc2;
if (Sw + eps < 1.0) {
pc2 = BrooksCorey::twoPhaseSatPcnw(*this, Sw + eps);
delta += eps;
}
else
pc2 = BrooksCorey::twoPhaseSatPcnw(*this, Sw);
if (Sw - eps > 0.0) {
pc1 = BrooksCorey::twoPhaseSatPcnw(*this, Sw - eps);
delta += eps;
}
else
pc1 = BrooksCorey::twoPhaseSatPcnw(*this, Sw);
return (pc2 - pc1)/delta;
}
Scalar pcnwLowSw_;
Scalar pcnwLow_;
Scalar pcnwSlopeLow_;
Scalar pcnwHigh_;
Scalar pcnwSlopeHigh_;
};
} // namespace Opm

199
opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp
View File

@ -113,8 +113,10 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
}
/*!
@ -124,7 +126,9 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = Sw(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
@ -135,8 +139,10 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
}
/*!
@ -151,11 +157,17 @@ public:
*
* \copydetails VanGenuchten::pC()
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
{ return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{
// retrieve the low and the high threshold saturations for the
// unregularized capillary pressure curve from the parameters
@ -206,14 +218,19 @@ public:
\copydetails VanGenuchten::Sw()
*
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& pC =
FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pC)
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pC)
{
// retrieve the low and the high threshold saturations for the
// unregularized capillary pressure curve from the parameters
@ -223,14 +240,13 @@ public:
// calculate the saturation which corrosponds to the
// saturation in the non-regularized verision of van
// Genuchten's law
Scalar Sw;
if (pC <= 0) {
// invert straight line for Sw > 1.0
Scalar m1 = params.pcnwSlopeHigh();
return pC/m1 + 1.0;
}
else
Sw = VanGenuchten::twoPhaseSatSw(params, pC);
Evaluation Sw = VanGenuchten::twoPhaseSatSw(params, pC);
// invert the regularization if necessary
if (Sw <= SwThLow) {
@ -243,8 +259,8 @@ public:
// invert spline between threshold saturation and 1.0
const Spline<Scalar>& spline = params.pcnwHighSpline();
return spline.intersectInterval(/*x0=*/SwThHigh, /*x1=*/1.0,
/*a=*/0.0, /*b=*/0.0, /*c=*/0.0, /*d=*/pC);
return spline.template intersectInterval<Evaluation>(/*x0=*/SwThHigh, /*x1=*/1.0,
/*a=*/0, /*b=*/0, /*c=*/0, /*d=*/pC);
}
return Sw;
@ -254,92 +270,14 @@ public:
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw<FluidState, Evaluation>(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pc)
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pc)
{ return 1 - twoPhaseSatSw(params, pc); }
/*!
* \brief A regularized version of the partial derivative
* of the \f$p_c(\overline S_w)\f$ w.r.t. effective saturation
* according to van Genuchten.
*
* regularized part:
* - low saturation: use the slope of the regularization point (i.e. no kink).
* - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line and use that slope (yes, there is a kink :-( ).
*
* For not-regularized part:
*
* \copydetails VanGenuchten::dpC_dSw()
*
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
// derivative of the regualarization
if (Sw < params.pcnwLowSw()) {
// the slope of the straight line used in pC()
return params.pcnwSlopeLow();
}
else if (params.pcnwHighSw() <= Sw) {
if (Sw < 1)
return params.pcnwHighSpline().evalDerivative(Sw);
else
// the slope of the straight line used for the right
// side of the capillary pressure function
return params.pcnwSlopeHigh();
}
return VanGenuchten::twoPhaseSatDPcnw_dSw(params, Sw);
}
/*!
* \brief A regularized version of the partial derivative
* of the \f$\overline S_w(p_c)\f$ w.r.t. cap.pressure
* according to van Genuchten.
*
* regularized part:
* - low saturation: use the slope of the regularization point (i.e. no kink).
* - high saturation: connect the high regularization point with \f$ \overline S_w =1\f$ by a straight line and use that slope (yes, there is a kink :-( ).
*
* For not-regularized part:
* \copydetails VanGenuchten::dSw_dpC()
*/
template <class FluidState>
static Scalar dSw_dpC(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
// calculate the saturation which corrosponds to the
// saturation in the non-regularized verision of van
// Genuchten's law
Scalar Sw;
if (pC < 0)
Sw = 1.5; // make sure we regularize below
else
Sw = VanGenuchten::Sw_raw(params, pC);
// derivative of the regularization
if (Sw < params.pcnwLowSw()) {
// same as in dpC_dSw() but inverted
return 1/params.pcnwSlopeLow();
}
if (Sw > params.pcnwHighSw()) {
if (Sw < 1)
return 1/params.pcnwHighSpline().evalDerivative(Sw);
// same as in dpC_dSw() but inverted
return 1/params.pcnwSlopHigh();
}
return VanGenuchten::dSw_dpnw(params, fs);
}
/*!
* \brief Regularized version of the relative permeability
* for the wetting phase of
@ -354,27 +292,27 @@ public:
* For not-regularized part:
\copydetails VanGenuchten::krw()
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
// regularize
if (Sw < 0)
return 0;
else if (Sw > 1)
return 1;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
return VanGenuchten::twoPhaseSatKrw(params, Sw);
const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
if (Sw <= 0.0 || Sw >= 1.0)
return 0.0;
typedef MathToolbox<Evaluation> Toolbox;
return VanGenuchten::twoPhaseSatDKrw_dSw(params, Sw);
// regularize
if (Sw <= 0)
return Toolbox::createConstant(0);
else if (Sw >= 1)
return Toolbox::createConstant(1);
return VanGenuchten::twoPhaseSatKrw(params, Sw);
}
/*!
@ -391,27 +329,28 @@ public:
\copydetails VanGenuchten::krn()
*
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
// regularize
if (Sw <= 0)
return 1;
else if (Sw >= 1)
return 0;
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
return VanGenuchten::twoPhaseSatKrn(params, Sw);
const Evaluation& Sw =
1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
if (Sw <= 0.0 || Sw >= 1.0)
return 0.0;
typedef MathToolbox<Evaluation> Toolbox;
return VanGenuchten::twoPhaseSatDKrn_dSw(params, Sw);
// regularize
if (Sw <= 0)
return Toolbox::createConstant(1);
else if (Sw >= 1)
return Toolbox::createConstant(0);
return VanGenuchten::twoPhaseSatKrn(params, Sw);
}
};

14
opm/material/fluidmatrixinteractions/RegularizedVanGenuchtenParams.hpp
View File

@ -71,11 +71,11 @@ public:
Parent::finalize();
pcnwLow_ = VanGenuchten::twoPhaseSatPcnw(*this, pcnwLowSw_);
pcnwSlopeLow_ = VanGenuchten::twoPhaseSatDPcnw_dSw(*this, pcnwLowSw_);
pcnwSlopeLow_ = dPcnw_dSw_(pcnwLowSw_);
pcnwHigh_ = VanGenuchten::twoPhaseSatPcnw(*this, pcnwHighSw_);
pcnwSlopeHigh_ = 2*(0.0 - pcnwHigh_)/(1.0 - pcnwHighSw_);
Scalar mThreshold = VanGenuchten::twoPhaseSatDPcnw_dSw(*this, pcnwHighSw_);
Scalar mThreshold = dPcnw_dSw_(pcnwHighSw_);
pcnwHighSpline_.set(pcnwHighSw_, 1.0, // x0, x1
pcnwHigh_, 0, // y0, y1
@ -164,6 +164,16 @@ private:
{ }
#endif
Scalar dPcnw_dSw_(Scalar Sw) const
{
// use finite differences to calculate the derivative w.r.t. Sw of the
// unregularized curve's capillary pressure.
const Scalar eps = 1e-7;
Scalar pc1 = VanGenuchten::twoPhaseSatPcnw(*this, Sw - eps);
Scalar pc2 = VanGenuchten::twoPhaseSatPcnw(*this, Sw + eps);
return (pc2 - pc1)/(2*eps);
}
Scalar pcnwLowSw_;
Scalar pcnwHighSw_;

118
opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp
View File

@ -88,10 +88,12 @@ public:
* \brief The capillary pressure-saturation curve.
*/
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
static void capillaryPressures(Container &values, const Params &params, const FluidState &fluidState)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fluidState);
}
/*!
@ -99,26 +101,31 @@ public:
* pressure differences.
*/
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
static void saturations(Container &values, const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not implemented: saturations()"); }
/*!
* \brief The relative permeabilities
*/
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fluidState)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
}
/*!
* \brief The capillary pressure-saturation curve
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fluidState)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatPcnw(params, Sw);
}
@ -126,89 +133,78 @@ public:
/*!
* \brief The saturation-capillary pressure curve
*/
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{ return params.pcnwSpline().eval(Sw, /*extrapolate=*/true); }
/*!
* \brief The saturation-capillary pressure curve
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
static Scalar twoPhaseSatSw(const Params &params, Scalar pC)
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pC)
{ OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); }
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fluidState)
{ return 1 - Sw<FluidState, Evaluation>(params, fluidState); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pC)
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pC)
{ return 1 - twoPhaseSatSw(params, pC); }
/*!
* \brief The partial derivative of the capillary pressure with
* regard to the saturation
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
assert(0 < Sw && Sw < 1);
return params.pcnwSpline().evalDerivative(Sw, /*extrapolate=*/true);
}
/*!
* \brief The relative permeability for the wetting phase of the
* porous medium
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, params.krwSpline().eval(Sw, /*extrapolate=*/true))); }
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(Traits::wettingPhaseIdx));
/*!
* \brief The derivative of the relative permeability of the
* wetting phase in regard to the wetting saturation of the
* porous medium
*/
template <class FluidState>
static Scalar dKrw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
return twoPhaseSatKrw(params, Sw);
}
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
{ return params.krwSpline().evalDerivative(Sw, /*extrapolate=*/true); }
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, params.krwSpline().eval(Sw, /*extrapolate=*/true)));
}
/*!
* \brief The relative permeability for the non-wetting phase
* of the porous medium
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
{ return std::max(0.0, std::min(1.0, params.krnSpline().eval(Sw, /*extrapolate=*/true))); }
const Evaluation& Sn =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(Traits::nonWettingPhaseIdx));
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the porous medium
*/
template <class FluidState>
static Scalar dKrn_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
return twoPhaseSatKrn(params, 1.0 - Sn);
}
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
{ return params.krnSpline().evalDerivative(Sw, /*extrapolate=*/true); }
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::max(0.0, Toolbox::min(1.0, params.krnSpline().eval(Sw, /*extrapolate=*/true)));
}
};
} // namespace Opm

153
opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp
View File

@ -27,6 +27,8 @@
#include "ThreePhaseParkerVanGenuchtenParams.hpp"
#include <opm/material/common/MathToolbox.hpp>
#include <algorithm>
namespace Opm {
@ -99,9 +101,11 @@ public:
const Params &params,
const FluidState &fluidState)
{
values[gasPhaseIdx] = pcgn(params, fluidState);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, fluidState);
values[nonWettingPhaseIdx] = 0;
values[wettingPhaseIdx] = - pcnw(params, fluidState);
values[wettingPhaseIdx] = - pcnw<FluidState, Evaluation>(params, fluidState);
}
/*!
@ -113,18 +117,20 @@ public:
* p_{c,gn} = p_g - p_n
* \f]
*/
template <class FluidState>
static Scalar pcgn(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcgn(const Params &params, const FluidState &fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
Scalar PC_VG_REG = 0.01;
// sum of liquid saturations
Scalar St =
fluidState.saturation(wettingPhaseIdx)
+ fluidState.saturation(nonWettingPhaseIdx);
const auto& St =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx))
+ FsToolbox::template toLhs<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
Scalar Se = (St - params.Swrx())/(1. - params.Swrx());
Evaluation Se = (St - params.Swrx())/(1. - params.Swrx());
// regularization
if (Se < 0.0)
@ -134,8 +140,8 @@ public:
if (Se>PC_VG_REG && Se<1-PC_VG_REG)
{
Scalar x = std::pow(Se,-1/params.vgM()) - 1;
return std::pow(x, 1 - params.vgM())/params.vgAlpha();
const Evaluation& x = Toolbox::pow(Se,-1/params.vgM()) - 1;
return Toolbox::pow(x, 1.0 - params.vgM())/params.vgAlpha();
}
// value and derivative at regularization point
@ -144,10 +150,10 @@ public:
Se_regu = PC_VG_REG;
else
Se_regu = 1-PC_VG_REG;
Scalar x = std::pow(Se_regu,-1/params.vgM())-1;
Scalar pc = std::pow(x, 1/params.vgN())/params.vgAlpha();
Scalar pc_prime =
std::pow(x, 1/params.vgN()-1)
const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
const Evaluation& pc = Toolbox::pow(x, 1.0/params.vgN())/params.vgAlpha();
const Evaluation& pc_prime =
Toolbox::pow(x, 1/params.vgN()-1)
* std::pow(Se_regu,-1/params.vgM()-1)
/ (-params.vgM())
/ params.vgAlpha()
@ -167,12 +173,15 @@ public:
* p_{c,nw} = p_n - p_w
* \f]
*/
template <class FluidState>
static Scalar pcnw(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fluidState)
{
Scalar Sw = fluidState.saturation(wettingPhaseIdx);
Scalar Se = (Sw-params.Swr())/(1.-params.Snr());
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx));
Evaluation Se = (Sw-params.Swr())/(1.-params.Snr());
Scalar PC_VG_REG = 0.01;
@ -183,8 +192,8 @@ public:
Se=1.0;
if (Se>PC_VG_REG && Se<1-PC_VG_REG) {
Scalar x = std::pow(Se,-1/params.vgM()) - 1.0;
x = std::pow(x, 1 - params.vgM());
Evaluation x = Toolbox::pow(Se,-1/params.vgM()) - 1.0;
x = Toolbox::pow(x, 1 - params.vgM());
return x/params.vgAlpha();
}
@ -195,10 +204,10 @@ public:
else
Se_regu = 1.0 - PC_VG_REG;
Scalar x = std::pow(Se_regu,-1/params.vgM())-1;
Scalar pc = std::pow(x, 1/params.vgN())/params.vgAlpha();
Scalar pc_prime =
std::pow(x,1/params.vgN()-1)
const Evaluation& x = std::pow(Se_regu,-1/params.vgM())-1;
const Evaluation& pc = Toolbox::pow(x, 1/params.vgN())/params.vgAlpha();
const Evaluation& pc_prime =
Toolbox::pow(x,1/params.vgN()-1)
* std::pow(Se_regu, -1.0/params.vgM() - 1)
/ (-params.vgM())
/ params.vgAlpha()
@ -222,25 +231,22 @@ public:
/*!
* \brief The saturation of the gas phase.
*/
template <class FluidState>
static Scalar Sg(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sg(const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not implemented: Sg()"); }
/*!
* \brief The saturation of the non-wetting (i.e., oil) phase.
*/
template <class FluidState>
static Scalar Sn(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not implemented: Sn()"); }
/*!
* \brief The saturation of the wetting (i.e., water) phase.
*/
template <class FluidState>
static Scalar Sw(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fluidState)
{ OPM_THROW(std::logic_error, "Not implemented: Sw()"); }
/*!
@ -251,9 +257,11 @@ public:
const Params &params,
const FluidState &fluidState)
{
values[wettingPhaseIdx] = krw(params, fluidState);
values[nonWettingPhaseIdx] = krn(params, fluidState);
values[gasPhaseIdx] = krg(params, fluidState);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
values[nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
}
/*!
@ -264,20 +272,23 @@ public:
* (see p61. in "Comparison of the Three-Phase Oil Relative Permeability Models"
* MOJDEH DELSHAD and GARY A. POPE, Transport in Porous Media 4 (1989), 59-83.)
*/
template <class FluidState>
static Scalar krw(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fluidState)
{
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx));
// transformation to effective saturation
Scalar Se = (fluidState.saturation(wettingPhaseIdx) - params.Swr()) / (1-params.Swr());
const Evaluation& Se = (Sw - params.Swr()) / (1-params.Swr());
// regularization
if(Se > 1.0) return 1.;
if(Se < 0.0) return 0.;
Scalar r = 1. - std::pow(1 - std::pow(Se, 1/params.vgM()), params.vgM());
return std::sqrt(Se)*r*r;
const Evaluation& r = 1. - Toolbox::pow(1 - Toolbox::pow(Se, 1/params.vgM()), params.vgM());
return Toolbox::sqrt(Se)*r*r;
}
/*!
@ -292,36 +303,39 @@ public:
* L. I. Oliveira, A. H. Demond, Journal of Contaminant Hydrology
* 66 (2003), 261-285
*/
template <class FluidState>
static Scalar krn(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fluidState)
{
Scalar Sn = fluidState.saturation(nonWettingPhaseIdx);
Scalar Sw = fluidState.saturation(wettingPhaseIdx);
Scalar Swe = std::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
Scalar Ste = std::min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sn =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(nonWettingPhaseIdx));
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(wettingPhaseIdx));
Evaluation Swe = Toolbox::min((Sw - params.Swr()) / (1 - params.Swr()), 1.);
Evaluation Ste = Toolbox::min((Sw + Sn - params.Swr()) / (1 - params.Swr()), 1.);
// regularization
if(Swe <= 0.0) Swe = 0.;
if(Ste <= 0.0) Ste = 0.;
if(Ste - Swe <= 0.0) return 0.;
Scalar krn_;
krn_ = std::pow(1 - std::pow(Swe, 1/params.vgM()), params.vgM());
krn_ -= std::pow(1 - std::pow(Ste, 1/params.vgM()), params.vgM());
Evaluation krn_;
krn_ = Toolbox::pow(1 - Toolbox::pow(Swe, 1/params.vgM()), params.vgM());
krn_ -= Toolbox::pow(1 - Toolbox::pow(Ste, 1/params.vgM()), params.vgM());
krn_ *= krn_;
if (params.krRegardsSnr())
{
// regard Snr in the permeability of the non-wetting
// phase, see Helmig1997
Scalar resIncluded =
std::max(std::min(Sw - params.Snr() / (1-params.Swr()), 1.0),
0.0);
krn_ *= std::sqrt(resIncluded );
const Evaluation& resIncluded =
Toolbox::max(Toolbox::min(Sw - params.Snr() / (1-params.Swr()), 1.0), 0.0);
krn_ *= Toolbox::sqrt(resIncluded );
}
else
krn_ *= std::sqrt(Sn / (1 - params.Swr()));
krn_ *= Toolbox::sqrt(Sn / (1 - params.Swr()));
return krn_;
}
@ -337,12 +351,15 @@ public:
* Three-Phase Oil Relative Permeability Models" M. Delshad and
* G. A. Pope, Transport in Porous Media 4 (1989), 59-83.)
*/
template <class FluidState>
static Scalar krg(const Params &params,
const FluidState &fluidState)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krg(const Params &params, const FluidState &fluidState)
{
Scalar Sg = fluidState.saturation(gasPhaseIdx);
Scalar Se = std::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
typedef MathToolbox<Evaluation> Toolbox;
const Evaluation& Sg =
FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));
const Evaluation& Se = Toolbox::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.);
// regularization
if(Se > 1.0)
@ -350,16 +367,16 @@ public:
if(Se < 0.0)
return 1.0;
Scalar scaleFactor = 1.;
Evaluation scaleFactor = 1.;
if (Sg<=0.1) {
scaleFactor = (Sg - params.Sgr())/(0.1 - params.Sgr());
if (scaleFactor < 0.)
scaleFactor = 0.;
return 0.0;
}
return scaleFactor
* std::pow(1 - Se, 1.0/3.)
* std::pow(1 - std::pow(Se, 1/params.vgM()), 2*params.vgM());
* Toolbox::pow(1 - Se, 1.0/3.)
* Toolbox::pow(1 - Toolbox::pow(Se, 1/params.vgM()), 2*params.vgM());
}
};
} // namespace Opm

174
opm/material/fluidmatrixinteractions/VanGenuchten.hpp
View File

@ -26,6 +26,8 @@
#include "VanGenuchtenParams.hpp"
#include <opm/material/common/MathToolbox.hpp>
#include <algorithm>
#include <cmath>
#include <cassert>
@ -108,8 +110,10 @@ public:
template <class Container, class FluidState>
static void capillaryPressures(Container &values, const Params &params, const FluidState &fs)
{
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = 0.0; // reference phase
values[Traits::nonWettingPhaseIdx] = pcnw(params, fs);
values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
}
/*!
@ -119,7 +123,9 @@ public:
template <class Container, class FluidState>
static void saturations(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = Sw(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
}
@ -136,8 +142,10 @@ public:
template <class Container, class FluidState>
static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
{
values[Traits::wettingPhaseIdx] = krw(params, fs);
values[Traits::nonWettingPhaseIdx] = krn(params, fs);
typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;
values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
}
/*!
@ -154,10 +162,14 @@ public:
* \param fs The fluid state for which the capillary pressure
* ought to be calculated
*/
template <class FluidState>
static Scalar pcnw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation pcnw(const Params &params, const FluidState &fs)
{
Scalar Sw = fs.saturation(Traits::wettingPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
assert(0 <= Sw && Sw <= 1);
return twoPhaseSatPcnw(params, Sw);
@ -177,8 +189,13 @@ public:
* required by the van Genuchten law.
* \param Sw The effective wetting phase saturation
*/
static Scalar twoPhaseSatPcnw(const Params &params, Scalar Sw)
{ return std::pow(std::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha(); }
template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
return Toolbox::pow(Toolbox::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha();
}
/*!
* \brief The saturation-capillary pressure curve according to van Genuchten.
@ -192,59 +209,39 @@ public:
* required by the van Genuchten law.
* \param fs The fluid state containing valid phase pressures
*/
template <class FluidState>
static Scalar Sw(const Params &params, const FluidState &fs)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sw(const Params &params, const FluidState &fs)
{
Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Evaluation pC =
FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
- FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
return twoPhaseSatSw(params, pC);
}
static Scalar twoPhaseSatSw(const Params &params, Scalar pC)
template <class Evaluation>
static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pC)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(pC >= 0);
return std::pow(std::pow(params.vgAlpha()*pC, params.vgN()) + 1, -params.vgM());
return Toolbox::pow(Toolbox::pow(params.vgAlpha()*pC, params.vgN()) + 1, -params.vgM());
}
/*!
* \brief Calculate the non-wetting phase saturations depending on
* the phase pressures.
*/
template <class FluidState>
static Scalar Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw(params, fs); }
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation Sn(const Params &params, const FluidState &fs)
{ return 1 - Sw<FluidState, Evaluation>(params, fs); }
static Scalar twoPhaseSatSn(const Params &params, Scalar pC)
template <class Evaluation>
static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pC)
{ return 1 - twoPhaseSatSw(params, pC); }
/*!
* \brief The partial derivative of the capillary pressure with
* regard to the saturation according to van Genuchten.
*
* This is equivalent to
* \f[
* \frac{\partial p_C}{\partial S_w} =
* -\frac{1}{\alpha n} ({S_w}^{-1/m} - 1)^{1/n - 1} {S_w}^{-1/m - 1}
* \f]
*
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state containing valid saturations
*/
template <class FluidState>
static Scalar dPcnw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDPcnw_dSw(const Params &params, Scalar Sw)
{
assert(0 < Sw && Sw < 1);
Scalar powSw = std::pow(Sw, -1/params.vgM());
return
- 1/(params.vgAlpha() * params.vgN() * Sw)
* std::pow(powSw - 1, 1/params.vgN() - 1) * powSw;
}
/*!
* \brief The relative permeability for the wetting phase of the
* medium according to van Genuchten's curve with Mualem
@ -255,43 +252,28 @@ public:
* \param fs The fluid state for which the relative permeability
* ought to be calculated
*/
template <class FluidState>
static Scalar krw(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatKrw(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krw(const Params &params, const FluidState &fs)
{
assert(0 <= Sw && Sw <= 1);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
Scalar r = 1. - std::pow(1 - std::pow(Sw, 1/params.vgM()), params.vgM());
return std::sqrt(Sw)*r*r;
const Evaluation& Sw =
FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
return twoPhaseSatKrw(params, Sw);
}
/*!
* \brief The derivative of the relative permeability of the
* wetting phase in regard to the wetting saturation of the
* medium implied according to the van Genuchten curve with
* Mualem parameters.
*
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the derivative
* ought to be calculated
*/
template <class FluidState>
static Scalar dKrw_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrw_dSw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
{
assert(0 <= Sw && Sw <= 1);
typedef MathToolbox<Evaluation> Toolbox;
const Scalar x = 1 - std::pow(Sw, 1.0/params.vgM());
const Scalar xToM = std::pow(x, params.vgM());
return (1 - xToM)/std::sqrt(Sw) * ( (1 - xToM)/2 + 2*xToM*(1-x)/x );
assert(0.0 <= Sw && Sw <= 1.0);
Evaluation r = 1.0 - Toolbox::pow(1.0 - Toolbox::pow(Sw, 1/params.vgM()), params.vgM());
return Toolbox::sqrt(Sw)*r*r;
}
/*!
* \brief The relative permeability for the non-wetting phase
* of the medium according to van Genuchten.
@ -301,43 +283,27 @@ public:
* \param fs The fluid state for which the derivative
* ought to be calculated
*/
template <class FluidState>
static Scalar krn(const Params &params, const FluidState &fs)
{ return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); }
static Scalar twoPhaseSatKrn(const Params &params, Scalar Sw)
template <class FluidState, class Evaluation = typename FluidState::Scalar>
static Evaluation krn(const Params &params, const FluidState &fs)
{
assert(0 <= Sw && Sw <= 1);
typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
return
std::pow(1 - Sw, 1.0/3) *
std::pow(1 - std::pow(Sw, 1/params.vgM()), 2*params.vgM());
const Evaluation& Sw =
1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
return twoPhaseSatKrn(params, Sw);
}
/*!
* \brief The derivative of the relative permeability for the
* non-wetting phase in regard to the wetting saturation of
* the medium as implied by the van Genuchten
* parameterization.
*
* \param params The parameter object expressing the coefficients
* required by the van Genuchten law.
* \param fs The fluid state for which the derivative
* ought to be calculated
*/
template <class FluidState>
static Scalar dKrn_dSw(const Params &params, const FluidState &fs)
{ return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); }
static Scalar twoPhaseSatDKrn_dSw(const Params &params, Scalar Sw)
template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params &params, Evaluation Sw)
{
typedef MathToolbox<Evaluation> Toolbox;
assert(0 <= Sw && Sw <= 1);
const Scalar x = std::pow(Sw, 1.0/params.vgM());
return
-std::pow(1 - x, 2*params.vgM())
*std::pow(1 - Sw, -2/3)
*(1.0/3 + 2*x/Sw);
Toolbox::pow(1 - Sw, 1.0/3) *
Toolbox::pow(1 - Toolbox::pow(Sw, 1/params.vgM()), 2*params.vgM());
}
};
} // namespace Opm

77
tests/test_fluidmatrixinteractions.cpp
View File

@ -25,6 +25,9 @@
*/
#include "config.h"
// include the local AD framwork
#include <opm/material/localad/Math.hpp>
// include all capillary pressure laws
#include <opm/material/fluidmatrixinteractions/BrooksCorey.hpp>
#include <opm/material/fluidmatrixinteractions/ParkerLenhard.hpp>
@ -122,10 +125,20 @@ void testGenericApi()
// test the generic methods which need to be implemented by
// all material laws
const FluidState fs;
double destValues[numPhases];
MaterialLaw::capillaryPressures(destValues, paramsConst, fs);
MaterialLaw::saturations(destValues, paramsConst, fs);
MaterialLaw::relativePermeabilities(destValues, paramsConst, fs);
{
double destValues[numPhases];
MaterialLaw::capillaryPressures(destValues, paramsConst, fs);
MaterialLaw::saturations(destValues, paramsConst, fs);
MaterialLaw::relativePermeabilities(destValues, paramsConst, fs);
}
{
typename FluidState::Scalar destValuesEval[numPhases];
MaterialLaw::capillaryPressures(destValuesEval, paramsConst, fs);
MaterialLaw::saturations(destValuesEval, paramsConst, fs);
MaterialLaw::relativePermeabilities(destValuesEval, paramsConst, fs);
}
}
}
@ -155,11 +168,19 @@ void testTwoPhaseApi()
const typename MaterialLaw::Params params;
OPM_UNUSED Scalar v;
v = MaterialLaw::pcnw(params, fs);
v = MaterialLaw::Sw(params, fs);
v = MaterialLaw::Sn(params, fs);
v = MaterialLaw::krw(params, fs);
v = MaterialLaw::krn(params, fs);
v = MaterialLaw::template pcnw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sn<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krn<FluidState, Scalar>(params, fs);
OPM_UNUSED typename FluidState::Scalar vEval;
vEval = MaterialLaw::pcnw(params, fs);
vEval = MaterialLaw::Sw(params, fs);
vEval = MaterialLaw::Sn(params, fs);
vEval = MaterialLaw::krw(params, fs);
vEval = MaterialLaw::krn(params, fs);
}
}
@ -194,6 +215,14 @@ void testTwoPhaseSatApi()
v = MaterialLaw::twoPhaseSatSn(params, Sw);
v = MaterialLaw::twoPhaseSatKrw(params, Sw);
v = MaterialLaw::twoPhaseSatKrn(params, Sw);
typename FluidState::Scalar SwEval = 0;
OPM_UNUSED typename FluidState::Scalar vEval;
vEval = MaterialLaw::twoPhaseSatPcnw(params, SwEval);
vEval = MaterialLaw::twoPhaseSatSw(params, SwEval);
vEval = MaterialLaw::twoPhaseSatSn(params, SwEval);
vEval = MaterialLaw::twoPhaseSatKrw(params, SwEval);
vEval = MaterialLaw::twoPhaseSatKrn(params, SwEval);
}
}
@ -217,13 +246,22 @@ void testThreePhaseApi()
const typename MaterialLaw::Params params;
OPM_UNUSED Scalar v;
v = MaterialLaw::pcnw(params, fs);
v = MaterialLaw::Sw(params, fs);
v = MaterialLaw::Sn(params, fs);
v = MaterialLaw::Sg(params, fs);
v = MaterialLaw::krw(params, fs);
v = MaterialLaw::krn(params, fs);
v = MaterialLaw::krg(params, fs);
v = MaterialLaw::template pcnw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sn<FluidState, Scalar>(params, fs);
v = MaterialLaw::template Sg<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krw<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krn<FluidState, Scalar>(params, fs);
v = MaterialLaw::template krg<FluidState, Scalar>(params, fs);
OPM_UNUSED typename FluidState::Scalar vEval;
vEval = MaterialLaw::pcnw(params, fs);
vEval = MaterialLaw::Sw(params, fs);
vEval = MaterialLaw::Sn(params, fs);
vEval = MaterialLaw::Sg(params, fs);
vEval = MaterialLaw::krw(params, fs);
vEval = MaterialLaw::krn(params, fs);
vEval = MaterialLaw::krg(params, fs);
}
}
@ -232,6 +270,8 @@ void testThreePhaseSatApi()
{
}
class TestAdTag;
int main(int argc, char **argv)
{
typedef double Scalar;
@ -253,8 +293,9 @@ int main(int argc, char **argv)
ThreePFluidSystem::oilPhaseIdx,
ThreePFluidSystem::gasPhaseIdx> ThreePhaseTraits;
typedef Opm::ImmiscibleFluidState<Scalar, TwoPFluidSystem> TwoPhaseFluidState;
typedef Opm::ImmiscibleFluidState<Scalar, ThreePFluidSystem> ThreePhaseFluidState;
typedef Opm::LocalAd::Evaluation<Scalar, TestAdTag, 3> Evaluation;
typedef Opm::ImmiscibleFluidState<Evaluation, TwoPFluidSystem> TwoPhaseFluidState;
typedef Opm::ImmiscibleFluidState<Evaluation, ThreePFluidSystem> ThreePhaseFluidState;
MyMpiHelper mpiHelper(argc, argv);