Mark Objects 'const' Where Possible

Also, don't use lifetime extension for 'Evaluation' objects where
these are constructed inside functions.  Doing so inhibits NRVO.
This commit is contained in:
Bård Skaflestad 2021-01-11 18:01:32 +01:00
parent f8f1bfc341
commit 247ce0e3fa
4 changed files with 67 additions and 58 deletions

View File

@ -220,7 +220,7 @@ public:
static Evaluation pcgn(const Params& params, static Evaluation pcgn(const Params& params,
const FluidState& fs) const FluidState& fs)
{ {
const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx)); const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw); return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
} }
@ -237,7 +237,7 @@ public:
static Evaluation pcnw(const Params& params, static Evaluation pcnw(const Params& params,
const FluidState& fs) const FluidState& fs)
{ {
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx)); const auto Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw); return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
} }
@ -316,7 +316,7 @@ public:
static Evaluation krg(const Params& params, static Evaluation krg(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx)); const Evaluation Sw = 1.0 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw); return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
} }
@ -327,7 +327,7 @@ public:
static Evaluation krw(const Params& params, static Evaluation krw(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)); const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw); return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
} }
@ -338,27 +338,28 @@ public:
static Evaluation krn(const Params& params, static Evaluation krn(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
Scalar Swco = params.Swl(); const Scalar Swco = params.Swl();
Evaluation Sw = const Evaluation Sw =
Opm::max(Evaluation(Swco), Opm::max(Evaluation(Swco),
Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx))); Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Evaluation Sw_ow = Sg + Sw; const Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Evaluation So_go = 1.0 - Sw_ow;
const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow); const Evaluation Sw_ow = Sg + Sw;
const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go); const Evaluation So_go = 1.0 - Sw_ow;
const Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
const Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
// avoid the division by zero: chose a regularized kro which is used if Sw - Swco // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
// < epsilon/2 and interpolate between the oridinary and the regularized kro between // < epsilon/2 and interpolate between the oridinary and the regularized kro between
// epsilon and epsilon/2 // epsilon and epsilon/2
const Scalar epsilon = 1e-5; const Scalar epsilon = 1e-5;
if (Opm::scalarValue(Sw_ow) - Swco < epsilon) { if (Opm::scalarValue(Sw_ow) - Swco < epsilon) {
Evaluation kro2 = (kro_ow + kro_go)/2;; const Evaluation kro2 = (kro_ow + kro_go)/2;;
if (Opm::scalarValue(Sw_ow) - Swco > epsilon/2) { if (Opm::scalarValue(Sw_ow) - Swco > epsilon/2) {
Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco); const Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2); const Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
return kro2*alpha + kro1*(1 - alpha); return kro2*alpha + kro1*(1 - alpha);
} }
@ -397,12 +398,12 @@ public:
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg); params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg);
} }
else { else {
Scalar Swco = params.Swl(); const Scalar Swco = params.Swl();
Sw = std::min(Scalar(1.0), std::max(Scalar(0.0), Sw)); Sw = std::min(Scalar(1.0), std::max(Scalar(0.0), Sw));
Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg)); Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
Scalar Sw_ow = Sg + std::max(Swco, Sw); const Scalar Sw_ow = Sg + std::max(Swco, Sw);
Scalar So_go = 1 + Sw_ow; const Scalar So_go = 1 + Sw_ow;
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow); params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg); params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);

View File

@ -151,16 +151,16 @@ public:
template <class Evaluation> template <class Evaluation>
static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& SwScaled) static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& SwScaled)
{ {
const Evaluation& SwUnscaled = scaledToUnscaledSatPc(params, SwScaled); const Evaluation SwUnscaled = scaledToUnscaledSatPc(params, SwScaled);
const Evaluation& pcUnscaled = EffLaw::twoPhaseSatPcnw(params.effectiveLawParams(), SwUnscaled); const Evaluation pcUnscaled = EffLaw::twoPhaseSatPcnw(params.effectiveLawParams(), SwUnscaled);
return unscaledToScaledPcnw_(params, pcUnscaled); return unscaledToScaledPcnw_(params, pcUnscaled);
} }
template <class Evaluation> template <class Evaluation>
static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnwScaled) static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnwScaled)
{ {
Evaluation pcnwUnscaled = scaledToUnscaledPcnw_(params, pcnwScaled); const Evaluation pcnwUnscaled = scaledToUnscaledPcnw_(params, pcnwScaled);
Evaluation SwUnscaled = EffLaw::twoPhaseSatPcnwInv(params.effectiveLawParams(), pcnwUnscaled); const Evaluation SwUnscaled = EffLaw::twoPhaseSatPcnwInv(params.effectiveLawParams(), pcnwUnscaled);
return unscaledToScaledSatPc(params, SwUnscaled); return unscaledToScaledSatPc(params, SwUnscaled);
} }
@ -223,16 +223,16 @@ public:
template <class Evaluation> template <class Evaluation>
static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& SwScaled) static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& SwScaled)
{ {
const Evaluation& SwUnscaled = scaledToUnscaledSatKrw(params, SwScaled); const Evaluation SwUnscaled = scaledToUnscaledSatKrw(params, SwScaled);
const Evaluation& krwUnscaled = EffLaw::twoPhaseSatKrw(params.effectiveLawParams(), SwUnscaled); const Evaluation krwUnscaled = EffLaw::twoPhaseSatKrw(params.effectiveLawParams(), SwUnscaled);
return unscaledToScaledKrw_(SwScaled, params, krwUnscaled); return unscaledToScaledKrw_(SwScaled, params, krwUnscaled);
} }
template <class Evaluation> template <class Evaluation>
static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krwScaled) static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krwScaled)
{ {
Evaluation krwUnscaled = scaledToUnscaledKrw_(params, krwScaled); const Evaluation krwUnscaled = scaledToUnscaledKrw_(params, krwScaled);
Evaluation SwUnscaled = EffLaw::twoPhaseSatKrwInv(params.effectiveLawParams(), krwUnscaled); const Evaluation SwUnscaled = EffLaw::twoPhaseSatKrwInv(params.effectiveLawParams(), krwUnscaled);
return unscaledToScaledSatKrw(params, SwUnscaled); return unscaledToScaledSatKrw(params, SwUnscaled);
} }
@ -248,16 +248,16 @@ public:
template <class Evaluation> template <class Evaluation>
static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& SwScaled) static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& SwScaled)
{ {
const Evaluation& SwUnscaled = scaledToUnscaledSatKrn(params, SwScaled); const Evaluation SwUnscaled = scaledToUnscaledSatKrn(params, SwScaled);
const Evaluation& krnUnscaled = EffLaw::twoPhaseSatKrn(params.effectiveLawParams(), SwUnscaled); const Evaluation krnUnscaled = EffLaw::twoPhaseSatKrn(params.effectiveLawParams(), SwUnscaled);
return unscaledToScaledKrn_(params, krnUnscaled); return unscaledToScaledKrn_(params, krnUnscaled);
} }
template <class Evaluation> template <class Evaluation>
static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krnScaled) static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krnScaled)
{ {
Evaluation krnUnscaled = scaledToUnscaledKrn_(params, krnScaled); const Evaluation krnUnscaled = scaledToUnscaledKrn_(params, krnScaled);
Evaluation SwUnscaled = EffLaw::twoPhaseSatKrnInv(params.effectiveLawParams(), krnUnscaled); const Evaluation SwUnscaled = EffLaw::twoPhaseSatKrnInv(params.effectiveLawParams(), krnUnscaled);
return unscaledToScaledSatKrn(params, SwUnscaled); return unscaledToScaledSatKrn(params, SwUnscaled);
} }

View File

@ -34,6 +34,8 @@
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <stdexcept>
#include <type_traits>
namespace Opm { namespace Opm {
@ -218,7 +220,7 @@ public:
static Evaluation pcgn(const Params& params, static Evaluation pcgn(const Params& params,
const FluidState& fs) const FluidState& fs)
{ {
const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx)); const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw); return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
} }
@ -235,10 +237,12 @@ public:
static Evaluation pcnw(const Params& params, static Evaluation pcnw(const Params& params,
const FluidState& fs) const FluidState& fs)
{ {
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx)); const auto Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
Valgrind::CheckDefined(Sw); Valgrind::CheckDefined(Sw);
const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
const auto result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
Valgrind::CheckDefined(result); Valgrind::CheckDefined(result);
return result; return result;
} }
@ -317,7 +321,7 @@ public:
static Evaluation krg(const Params& params, static Evaluation krg(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx)); const Evaluation Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw); return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
} }
@ -328,7 +332,7 @@ public:
static Evaluation krw(const Params& params, static Evaluation krw(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)); const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw); return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
} }
@ -342,16 +346,16 @@ public:
// the Eclipse docu is inconsistent with naming the variable of connate water: In // the Eclipse docu is inconsistent with naming the variable of connate water: In
// some places the connate water saturation is represented by "Swl", in others // some places the connate water saturation is represented by "Swl", in others
// "Swco" is used. // "Swco" is used.
Scalar Swco = params.Swl(); const Scalar Swco = params.Swl();
// oil relperm at connate water saturations (with Sg=0) // oil relperm at connate water saturations (with Sg=0)
Scalar krocw = params.krocw(); const Scalar krocw = params.krocw();
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)); const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx)); const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw); const Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Swco); const Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Swco);
Evaluation beta; Evaluation beta;
if (Sw <= Swco) if (Sw <= Swco)
@ -360,9 +364,9 @@ public:
// there seems to be an error in the ECL documentation: using the approach to // there seems to be an error in the ECL documentation: using the approach to
// the scaled saturations as described there leads to significant deviations // the scaled saturations as described there leads to significant deviations
// from the results produced by Eclipse 100. // from the results produced by Eclipse 100.
Evaluation SSw = (Sw - Swco)/(1.0 - Swco); const Evaluation SSw = (Sw - Swco)/(1.0 - Swco);
Evaluation SSg = Sg/(1.0 - Swco); const Evaluation SSg = Sg/(1.0 - Swco);
Evaluation SSo = 1.0 - SSw - SSg; const Evaluation SSo = 1.0 - SSw - SSg;
if (SSw >= 1.0 || SSg >= 1.0) if (SSw >= 1.0 || SSg >= 1.0)
beta = 1.0; beta = 1.0;
@ -383,8 +387,8 @@ public:
template <class FluidState> template <class FluidState>
static void updateHysteresis(Params& params, const FluidState& fluidState) static void updateHysteresis(Params& params, const FluidState& fluidState)
{ {
Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx)); const Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx)); const Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw); params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg); params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg);

View File

@ -33,6 +33,8 @@
#include <opm/material/common/MathToolbox.hpp> #include <opm/material/common/MathToolbox.hpp>
#include <algorithm> #include <algorithm>
#include <stdexcept>
#include <type_traits>
namespace Opm { namespace Opm {
@ -219,7 +221,7 @@ public:
static Evaluation pcgn(const Params& params, static Evaluation pcgn(const Params& params,
const FluidState& fs) const FluidState& fs)
{ {
const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx)); const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw); return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
} }
@ -236,10 +238,12 @@ public:
static Evaluation pcnw(const Params& params, static Evaluation pcnw(const Params& params,
const FluidState& fs) const FluidState& fs)
{ {
const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx)); const auto Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
Valgrind::CheckDefined(Sw); Valgrind::CheckDefined(Sw);
const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
const auto result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
Valgrind::CheckDefined(result); Valgrind::CheckDefined(result);
return result; return result;
} }
@ -318,7 +322,7 @@ public:
static Evaluation krg(const Params& params, static Evaluation krg(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx)); const Evaluation Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw); return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
} }
@ -329,7 +333,7 @@ public:
static Evaluation krw(const Params& params, static Evaluation krw(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)); const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw); return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
} }
@ -340,15 +344,15 @@ public:
static Evaluation krn(const Params& params, static Evaluation krn(const Params& params,
const FluidState& fluidState) const FluidState& fluidState)
{ {
Scalar Swco = params.Swl(); const Scalar Swco = params.Swl();
const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)); const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx)); const Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco); const Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco);
Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw); const Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
Evaluation krw = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw); const Evaluation krw = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Sg); const Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Sg);
Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg); const Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg);
return krocw * ((krow/krocw + krw) * (krog/krocw + krg) - krw - krg); return krocw * ((krow/krocw + krw) * (krog/krocw + krg) - krw - krg);
} }