diff --git a/opm/material/fluidmatrixinteractions/BrooksCorey.hpp b/opm/material/fluidmatrixinteractions/BrooksCorey.hpp index bfaf44a45..21bdf1dcf 100644 --- a/opm/material/fluidmatrixinteractions/BrooksCorey.hpp +++ b/opm/material/fluidmatrixinteractions/BrooksCorey.hpp @@ -27,6 +27,8 @@ #include "BrooksCoreyParams.hpp" +#include + #include #include #include @@ -92,8 +94,10 @@ public: template static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) { + typedef typename std::remove_reference::type Evaluation; + values[Traits::wettingPhaseIdx] = 0.0; // reference phase - values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); + values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); } /*! @@ -103,7 +107,9 @@ public: template static void saturations(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = Sw(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = Sw(params, fs); values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx]; } @@ -120,8 +126,10 @@ public: template static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = krw(params, fs); - values[Traits::nonWettingPhaseIdx] = krn(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = krw(params, fs); + values[Traits::nonWettingPhaseIdx] = krn(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 - static Scalar pcnw(const Params ¶ms, const FluidState &fs) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) { - Scalar Sw = fs.saturation(Traits::wettingPhaseIdx); + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = + FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + + assert(0 <= Sw && Sw <= 1); + return twoPhaseSatPcnw(params, Sw); } - static Scalar twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) { + typedef MathToolbox 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 - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) { - Scalar pc = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx); - return twoPhaseSatSw(params, pc); + typedef MathToolbox FsToolbox; + + Evaluation pC = + FsToolbox::template toLhs(fs.pressure(Traits::nonWettingPhaseIdx)) + - FsToolbox::template toLhs(fs.pressure(Traits::wettingPhaseIdx)); + return twoPhaseSatSw(params, pC); } - static Scalar twoPhaseSatSw(const Params ¶ms, Scalar pc) + template + static Evaluation twoPhaseSatSw(const Params ¶ms, const Evaluation& pc) { + typedef MathToolbox 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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return 1 - Sw(params, fs); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) + { return 1 - Sw(params, fs); } - static Scalar twoPhaseSatSn(const Params ¶ms, Scalar pc) + template + static Evaluation twoPhaseSatSn(const Params ¶ms, 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 - static Scalar dPcnw_dSw(const Params ¶ms, const FluidState &fs) - { - Scalar Sw = fs.saturation(Traits::wettingPhaseIdx); - return twoPhaseSatDPcnw_dSw(params, Sw); - } - - static Scalar twoPhaseSatDPcnw_dSw(const Params ¶ms, 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 - static Scalar dSw_dpcnw(const Params ¶ms, const FluidState &fs) - { - Scalar Sw = fs.saturation(Traits::wettingPhaseIdx); - return twoPhaseSatDSw_dpcnw(params, Sw); - } - - static Scalar twoPhaseSatDSw_dpcnw(const Params ¶ms, 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 - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatKrw(const Params ¶ms, Scalar Sw) + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) { - assert(0 <= Sw && Sw <= 1); + typedef MathToolbox FsToolbox; - return std::pow(Sw, 2.0/params.lambda() + 3); + const auto& Sw = + FsToolbox::template toLhs(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 ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) { + typedef MathToolbox 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 - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); } - - static Scalar twoPhaseSatKrn(const Params ¶ms, Scalar Sw) + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) { + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = + 1.0 - FsToolbox::template toLhs(fs.saturation(Traits::nonWettingPhaseIdx)); + + return twoPhaseSatKrn(params, Sw); + } + + template + static Evaluation twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox 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 ¶ms, 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 diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp index e9485bb54..1ea17d6a0 100644 --- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp @@ -26,6 +26,7 @@ #include "EclDefaultMaterialParams.hpp" #include +#include #include #include @@ -133,9 +134,13 @@ public: const Params ¶ms, const FluidState &state) { - values[gasPhaseIdx] = pcgn(params, state); + typedef typename std::remove_reference::type Evaluation; + values[gasPhaseIdx] = pcgn(params, state); values[oilPhaseIdx] = 0; - values[waterPhaseIdx] = - pcnw(params, state); + values[waterPhaseIdx] = - pcnw(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 - static Scalar pcgn(const Params ¶ms, - const FluidState &fs) + template + static Evaluation pcgn(const Params ¶ms, + const FluidState &fs) { - Scalar Sw = 1 - fs.saturation(gasPhaseIdx); + typedef MathToolbox FsToolbox; + + const auto& Sw = 1.0 - FsToolbox::template toLhs(fs.saturation(gasPhaseIdx)); return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw); } @@ -164,12 +171,17 @@ public: * p_{c,nw} = p_n - p_w * \f] */ - template - static Scalar pcnw(const Params ¶ms, - const FluidState &fs) + template + static Evaluation pcnw(const Params ¶ms, + const FluidState &fs) { - Scalar Sw = fs.saturation(waterPhaseIdx); - return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw); + typedef MathToolbox FsToolbox; + + const auto& Sw = FsToolbox::template toLhs(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 - static Scalar Sg(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation Sg(const Params ¶ms, + 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 - static Scalar Sn(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation Sn(const Params ¶ms, + 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 - static Scalar Sw(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation Sw(const Params ¶ms, + const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not implemented: Sw()"); } @@ -233,57 +245,69 @@ public: const Params ¶ms, const FluidState &fluidState) { - values[waterPhaseIdx] = krw(params, fluidState); - values[oilPhaseIdx] = krn(params, fluidState); - values[gasPhaseIdx] = krg(params, fluidState); + typedef typename std::remove_reference::type Evaluation; + + values[waterPhaseIdx] = krw(params, fluidState); + values[oilPhaseIdx] = krn(params, fluidState); + values[gasPhaseIdx] = krg(params, fluidState); } /*! * \brief The relative permeability of the gas phase. */ - template - static Scalar krg(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation krg(const Params ¶ms, + const FluidState &fluidState) { - Scalar Sw = 1 - fluidState.saturation(gasPhaseIdx); + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = 1 - FsToolbox::template toLhs(fluidState.saturation(gasPhaseIdx)); return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw); } /*! * \brief The relative permeability of the wetting phase. */ - template - static Scalar krw(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation krw(const Params ¶ms, + const FluidState &fluidState) { - Scalar Sw = fluidState.saturation(waterPhaseIdx); + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = FsToolbox::template toLhs(fluidState.saturation(waterPhaseIdx)); return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw); } /*! * \brief The relative permeability of the non-wetting (i.e., oil) phase. */ - template - static Scalar krn(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation krn(const Params ¶ms, + const FluidState &fluidState) { + typedef MathToolbox Toolbox; + typedef MathToolbox 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(fluidState.saturation(waterPhaseIdx)); +// Evaluation So = FsToolbox::template toLhs(fluidState.saturation(oilPhaseIdx)); + Evaluation Sg = FsToolbox::template toLhs(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)); } } }; diff --git a/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp b/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp index ac2998209..464209bfd 100644 --- a/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp +++ b/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp @@ -122,7 +122,7 @@ public: phaseIdx)); } - EffLaw::capillaryPressures(values, params, overlayFs); + EffLaw::template capillaryPressures(values, params, overlayFs); } /*! @@ -148,7 +148,7 @@ public: phaseIdx)); } - EffLaw::relativePermeabilities(values, params, overlayFs); + EffLaw::template relativePermeabilities(values, params, overlayFs); } /*! @@ -162,8 +162,8 @@ public: * constitutive relation (e.g. Brooks & Corey, van * Genuchten, linear...) */ - template - static Scalar pcnw(const Params ¶ms, const FluidState &fs) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) { typedef Opm::SaturationOverlayFluidState OverlayFluidState; @@ -179,14 +179,14 @@ public: phaseIdx)); } - return EffLaw::pcnw(params, overlayFs); + return EffLaw::template pcnw(params, overlayFs); } - template - static typename std::enable_if::type - twoPhaseSatPcnw(const Params ¶ms, Scalar SwAbs) + template + static typename std::enable_if::type + twoPhaseSatPcnw(const Params ¶ms, 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 static void saturations(Container &values, const Params ¶ms, const FluidState &fs) { - EffLaw::saturations(values, params, fs); + EffLaw::template saturations(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 - static Scalar Sw(const Params ¶ms, const FluidState &fs) - { return absoluteSaturation(params, EffLaw::Sw(params, fs), Traits::wettingPhaseIdx); } + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) + { + return absoluteSaturation(params, + EffLaw::template Sw(params, fs), + Traits::wettingPhaseIdx); + } - template - static typename std::enable_if::type - twoPhaseSatSw(const Params ¶ms, Scalar Sw) - { return absoluteSaturation(params, EffLaw::twoPhaseSatSw(params, Sw), Traits::wettingPhaseIdx); } + template + static typename std::enable_if::type + twoPhaseSatSw(const Params ¶ms, 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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return absoluteSaturation(params, EffLaw::Sn(params, fs), Traits::nonWettingPhaseIdx); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) + { + return absoluteSaturation(params, + EffLaw::template Sn(params, fs), + Traits::nonWettingPhaseIdx); + } - template - static typename std::enable_if::type - twoPhaseSatSn(const Params ¶ms, Scalar Sw) - { return absoluteSaturation(params, EffLaw::twoPhaseSatSn(params, Sw), Traits::nonWettingPhaseIdx); } + template + static typename std::enable_if::type + twoPhaseSatSn(const Params ¶ms, 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 - static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type + template + static typename std::enable_if< (Traits::numPhases > 2), Evaluation>::type Sg(const Params ¶ms, 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 ¶ms, 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 ¶ms, Scalar pC) - { - return EffLaw::dSw_dpC(params, pC)*dSw_dSwe_(params); + return absoluteSaturation(params, + EffLaw::template Sg(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 - static Scalar krw(const Params ¶ms, const FluidState &fs) + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) { typedef Opm::SaturationOverlayFluidState OverlayFluidState; @@ -308,19 +284,19 @@ public: phaseIdx)); } - return EffLaw::krw(params, overlayFs); + return EffLaw::template krw(params, overlayFs); } - template - static typename std::enable_if::type - twoPhaseSatKrw(const Params ¶ms, Scalar Sw) + template + static typename std::enable_if::type + twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) { return EffLaw::twoPhaseSatKrw(params, effectiveSaturation(params, Sw, Traits::nonWettingPhaseIdx)); } /*! * \brief The relative permeability of the non-wetting phase. */ - template - static Scalar krn(const Params ¶ms, const FluidState &fs) + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) { typedef Opm::SaturationOverlayFluidState OverlayFluidState; @@ -336,12 +312,12 @@ public: phaseIdx)); } - return EffLaw::krn(params, overlayFs); + return EffLaw::template krn(params, overlayFs); } - template - static typename std::enable_if::type - twoPhaseSatKrn(const Params ¶ms, Scalar Sw) + template + static typename std::enable_if::type + twoPhaseSatKrn(const Params ¶ms, 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 - static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type + template + static typename std::enable_if< (Traits::numPhases > 2), Evaluation>::type krg(const Params ¶ms, const FluidState &fs) { typedef Opm::SaturationOverlayFluidState OverlayFluidState; @@ -367,20 +343,22 @@ public: phaseIdx)); } - return EffLaw::krg(params, overlayFs); + return EffLaw::template krg(params, overlayFs); } /*! * \brief Convert an absolute saturation to an effective one. */ - static Scalar effectiveSaturation(const Params ¶ms, Scalar S, int phaseIdx) - { return (S - params.residualSaturation(phaseIdx))/(1 - params.sumResidualSaturations()); } + template + static Evaluation effectiveSaturation(const Params ¶ms, 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 ¶ms, Scalar S, int phaseIdx) - { return S*(1 - params.sumResidualSaturations()) + params.residualSaturation(phaseIdx); } + template + static Evaluation absoluteSaturation(const Params ¶ms, const Evaluation& S, int phaseIdx) + { return S*(1.0 - params.sumResidualSaturations()) + params.residualSaturation(phaseIdx); } private: /*! diff --git a/opm/material/fluidmatrixinteractions/LinearMaterial.hpp b/opm/material/fluidmatrixinteractions/LinearMaterial.hpp index 61a913c82..f81f13e22 100644 --- a/opm/material/fluidmatrixinteractions/LinearMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/LinearMaterial.hpp @@ -25,8 +25,8 @@ #include "LinearMaterialParams.hpp" +#include #include - #include #include @@ -97,8 +97,12 @@ public: const Params ¶ms, const FluidState &state) { + typedef typename std::remove_reference::type Evaluation; + typedef MathToolbox FsToolbox; + for (int phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) { - Scalar S = state.saturation(phaseIdx); + const Evaluation& S = + FsToolbox::template toLhs(state.saturation(phaseIdx)); Valgrind::CheckDefined(S); values[phaseIdx] = @@ -126,47 +130,55 @@ public: const Params ¶ms, const FluidState &state) { + typedef typename std::remove_reference::type Evaluation; + typedef MathToolbox Toolbox; + typedef MathToolbox FsToolbox; + for (int phaseIdx = 0; phaseIdx < Traits::numPhases; ++phaseIdx) { - Scalar S = state.saturation(phaseIdx); + const Evaluation& S = + FsToolbox::template toLhs(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 - static Scalar pcnw(const Params ¶ms, const FluidState &fs) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) { - Scalar S = fs.saturation(Traits::wettingPhaseIdx); - Valgrind::CheckDefined(S); + typedef MathToolbox FsToolbox; + const Evaluation& Sw = + FsToolbox::template toLhs(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(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 - static typename std::enable_if::type - twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + template + static typename std::enable_if::type + twoPhaseSatPcnw(const Params ¶ms, 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 - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) { OPM_THROW(std::runtime_error, "Not implemented: Sw()"); } - template - static typename std::enable_if::type - twoPhaseSatSw(const Params ¶ms, Scalar Sw) + template + static typename std::enable_if::type + twoPhaseSatSw(const Params ¶ms, 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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) { OPM_THROW(std::runtime_error, "Not implemented: Sn()"); } - template - static typename std::enable_if::type - twoPhaseSatSn(const Params ¶ms, Scalar Sw) + template + static typename std::enable_if::type + twoPhaseSatSn(const Params ¶ms, 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 - static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type + template + static typename std::enable_if::type Sg(const Params ¶ms, const FluidState &fs) { OPM_THROW(std::runtime_error, "Not implemented: Sg()"); } /*! * \brief The relative permability of the wetting phase */ - template - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return std::max(0.0, std::min(1.0, fs.saturation(Traits::wettingPhaseIdx))); } + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) + { + typedef MathToolbox Toolbox; + typedef MathToolbox FsToolbox; - template - static typename std::enable_if::type - twoPhaseSatKrw(const Params ¶ms, Scalar Sw) - { return std::max(0.0, std::min(1.0, Sw)); } + const Evaluation& Sw = + FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + return Toolbox::max(0.0, Toolbox::min(1.0, Sw)); + } + + template + static typename std::enable_if::type + twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox Toolbox; + + return Toolbox::max(0.0, Toolbox::min(1.0, Sw)); + } /*! * \brief The relative permability of the liquid non-wetting phase */ - template - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return std::max(0.0, std::min(1.0, fs.saturation(Traits::nonWettingPhaseIdx))); } + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) + { + typedef MathToolbox Toolbox; + typedef MathToolbox FsToolbox; - template - static typename std::enable_if::type - twoPhaseSatKrn(const Params ¶ms, Scalar Sw) - { return std::max(0.0, std::min(1.0, 1 - Sw)); } + const Evaluation& Sn = + FsToolbox::template toLhs(fs.saturation(Traits::nonWettingPhaseIdx)); + return Toolbox::max(0.0, Toolbox::min(1.0, Sn)); + } + + template + static typename std::enable_if::type + twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox 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 - static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type + template + static typename std::enable_if::type krg(const Params ¶ms, const FluidState &fs) - { return std::max(0.0, std::min(1.0, fs.saturation(Traits::gasPhaseIdx))); } + { + typedef MathToolbox Toolbox; + typedef MathToolbox FsToolbox; + + const Evaluation& Sg = + FsToolbox::template toLhs(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 - static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type + template + static typename std::enable_if::type pcgn(const Params ¶ms, const FluidState &fs) { - Scalar S = fs.saturation(Traits::nonWettingPhaseIdx); - Valgrind::CheckDefined(S); + typedef MathToolbox FsToolbox; - Scalar nPhasePressure = - S*params.pcMaxSat(Traits::nonWettingPhaseIdx) + - (1.0 - S)*params.pcMinSat(Traits::nonWettingPhaseIdx); + const Evaluation& Sn = + FsToolbox::template toLhs(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(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; } diff --git a/opm/material/fluidmatrixinteractions/MaterialTraits.hpp b/opm/material/fluidmatrixinteractions/MaterialTraits.hpp index 0d5f73d87..030530348 100644 --- a/opm/material/fluidmatrixinteractions/MaterialTraits.hpp +++ b/opm/material/fluidmatrixinteractions/MaterialTraits.hpp @@ -35,13 +35,16 @@ namespace Opm { * * This traits class is intended to be used by the NullMaterial */ -template +template 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 +template 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 +template 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; diff --git a/opm/material/fluidmatrixinteractions/NullMaterial.hpp b/opm/material/fluidmatrixinteractions/NullMaterial.hpp index 3faaea6e0..78bb0ec70 100644 --- a/opm/material/fluidmatrixinteractions/NullMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/NullMaterial.hpp @@ -27,6 +27,7 @@ #include #include +#include #include @@ -107,50 +108,56 @@ public: template static void relativePermeabilities(ContainerT &values, const Params ¶ms, - const FluidState &state) + const FluidState &fluidState) { + typedef typename std::remove_reference::type Evaluation; + typedef MathToolbox FsToolbox; + typedef MathToolbox 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(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 - static typename std::enable_if<(numPhases > 1), ScalarT>::type - pcnw(const Params ¶ms, const FluidState &fs) + template + static typename std::enable_if<(numPhases > 1), Evaluation>::type + pcnw(const Params ¶ms, const FluidState &fluidState) { return 0; } - template - static typename std::enable_if::type - twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + template + static typename std::enable_if::type + twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) { return 0; } /*! * \brief Calculate wetting phase saturation given that the rest * of the fluid state has been initialized */ - template - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Scalar Sw(const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not defined: Sw()"); } - template - static typename std::enable_if::type - twoPhaseSatSw(const Params ¶ms, Scalar pcnw) + template + static typename std::enable_if::type + twoPhaseSatSw(const Params ¶ms, 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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) + template + static Scalar Sn(const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not defined: Sn()"); } - template - static typename std::enable_if::type - twoPhaseSatSn(const Params ¶ms, Scalar pcnw) + template + static typename std::enable_if::type + twoPhaseSatSn(const Params ¶ms, 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 - static typename std::enable_if< (numPhases > 2), ScalarT>::type - Sg(const Params ¶ms, const FluidState &fs) + template + static typename std::enable_if< (numPhases > 2), Evaluation>::type + Sg(const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not defined: Sg()"); } /*! * \brief The relative permability of the wetting phase */ - template - static typename std::enable_if<(numPhases > 1), ScalarT>::type - krw(const Params ¶ms, const FluidState &fs) - { return std::max(0.0, std::min(1.0, fs.saturation(Traits::wettingPhaseIdx))); } + template + static typename std::enable_if<(numPhases > 1), Evaluation>::type + krw(const Params ¶ms, const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; + typedef MathToolbox Toolbox; - template - static typename std::enable_if::type - twoPhaseSatKrw(const Params ¶ms, Scalar Sw) - { return std::max(0.0, std::min(1.0, Sw)); } + const Evaluation& Sw = + FsToolbox::template toLhs(fluidState.saturation(Traits::wettingPhaseIdx)); + + return Toolbox::max(0.0, Toolbox::min(1.0, Sw)); + } + + template + static typename std::enable_if::type + twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox Toolbox; + + return Toolbox::max(0.0, Toolbox::min(1.0, Sw)); + } /*! * \brief The relative permability of the liquid non-wetting phase */ - template - static typename std::enable_if<(numPhases > 1), ScalarT>::type - krn(const Params ¶ms, const FluidState &fs) - { return std::max(0.0, std::min(1.0, fs.saturation(Traits::nonWettingPhaseIdx))); } + template + static typename std::enable_if<(numPhases > 1), Evaluation>::type + krn(const Params ¶ms, const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; + typedef MathToolbox Toolbox; - template - static typename std::enable_if::type - twoPhaseSatKrn(const Params ¶ms, Scalar Sw) - { return std::max(0.0, std::min(1.0, 1 - Sw)); } + const Evaluation& Sn = + FsToolbox::template toLhs(fluidState.saturation(Traits::nonWettingPhaseIdx)); + + return Toolbox::max(0.0, Toolbox::min(1.0, Sn)); + } + + template + static typename std::enable_if::type + twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox 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 - static typename std::enable_if< (numPhases > 2), ScalarT>::type - krg(const Params ¶ms, const FluidState &fs) - { return std::max(0.0, std::min(1.0, fs.saturation(Traits::gasPhaseIdx))); } + template + static typename std::enable_if< (numPhases > 2), Evaluation>::type + krg(const Params ¶ms, const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; + typedef MathToolbox Toolbox; + + const Evaluation& Sg = + FsToolbox::template toLhs(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 - static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type - pcgn(const Params ¶ms, const FluidState &fs) + template + static typename std::enable_if< (Traits::numPhases > 2), Evaluation>::type + pcgn(const Params ¶ms, const FluidState &fluidState) { return 0; } }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/ParkerLenhard.hpp b/opm/material/fluidmatrixinteractions/ParkerLenhard.hpp index cd10ad78f..7b342e325 100644 --- a/opm/material/fluidmatrixinteractions/ParkerLenhard.hpp +++ b/opm/material/fluidmatrixinteractions/ParkerLenhard.hpp @@ -299,7 +299,9 @@ public: template static void update(Params ¶ms, const FluidState &fs) { - Scalar Sw = fs.saturation(Traits::wettingPhaseIdx); + typedef MathToolbox 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(params, fs); Scalar Sw_mic = VanGenuchten::twoPhaseSatSw(params.micParams(), pc); Scalar Sw_mdc = VanGenuchten::twoPhaseSatSw(params.mdcParams(), pc); @@ -339,8 +341,10 @@ public: template static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) { + typedef typename std::remove_reference::type Evaluation; + values[Traits::wettingPhaseIdx] = 0.0; // reference phase - values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); + values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); } /*! @@ -358,25 +362,37 @@ public: template static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = krw(params, fs); - values[Traits::nonWettingPhaseIdx] = krn(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = krw(params, fs); + values[Traits::nonWettingPhaseIdx] = krn(params, fs); } /*! * \brief Returns the capillary pressure dependend on * the phase saturations. */ - template - static Scalar pcnw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) { + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = + FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + + return twoPhaseSatPcnw(params, Sw); + } + + template + static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox 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 ¶ms, Scalar Sw) - { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatDPcnw_dSw()"); } - /*! * \brief Calculate the wetting phase saturations depending on * the phase pressures. */ - template - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) { OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::Sw()"); } - static Scalar twoPhaseSatSw(const Params ¶ms, Scalar pc) + template + static Evaluation twoPhaseSatSw(const Params ¶ms, 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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return 1 - Sw(params, fs); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) + { return 1 - Sw(params, fs); } - static Scalar twoPhaseSatSn(const Params ¶ms, Scalar pc) + template + static Evaluation twoPhaseSatSn(const Params ¶ms, const Evaluation& pc) { OPM_THROW(std::logic_error, "Not implemented: ParkerLenhard::twoPhaseSatSn()"); } /*! * \brief The relative permeability for the wetting phase of * the medium. */ - template - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); } + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) + { + typedef MathToolbox FsToolbox; - static Scalar twoPhaseSatKrw(const Params ¶ms, Scalar Sw) + const Evaluation& Sw = + FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + + return twoPhaseSatKrw(params, Sw); + } + + template + static Evaluation twoPhaseSatKrw(const Params ¶ms, 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 ¶ms, Scalar Sw) - { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatDKrw_dSw()"); } - /*! * \brief The relative permeability for the non-wetting phase * of the params. */ - template - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrn(params, fs.saturation(Traits::wettingPhaseIdx)); } + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) + { + typedef MathToolbox FsToolbox; - static Scalar twoPhaseSatKrn(const Params ¶ms, Scalar Sw) + const Evaluation& Sw = + FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + + return twoPhaseSatKrn(params, Sw); + } + + template + static Evaluation twoPhaseSatKrn(const Params ¶ms, 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 ¶ms, 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 ¶ms, Scalar Sw) + template + static Evaluation absoluteToApparentSw_(const Params ¶ms, 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 ¶ms, Scalar Sw) + template + static Evaluation absoluteToEffectiveSw_(const Params ¶ms, 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 ¶ms, Scalar Swe) + template + static Evaluation effectiveToAbsoluteSw_(const Params ¶ms, 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 ¶ms, Scalar Sw) + template + static Evaluation computeCurrentSnr_(const Params ¶ms, 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 ¶ms, Scalar Sw) + template + static Evaluation trappedEffectiveSn_(const Params ¶ms, 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 ¶ms, Scalar Swe) + template + static Evaluation effectiveToApparentSw_(const Params ¶ms, 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 ¶ms, Scalar Swapp) + template + static Evaluation apparentToEffectiveSw_(const Params ¶ms, const Evaluation& Swapp) { Scalar SwePisc = absoluteToEffectiveSw_(params, params.pisc()->Sw()); if (params.pisc() == NULL || Swapp <= SwePisc) { diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp index 2e54dbdc7..d188f59fd 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterial.hpp @@ -27,6 +27,7 @@ #include #include +#include #include #include @@ -94,8 +95,10 @@ public: template static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) { + typedef typename std::remove_reference::type Evaluation; + values[Traits::wettingPhaseIdx] = 0.0; // reference phase - values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); + values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); } /*! @@ -112,17 +115,21 @@ public: template static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = krw(params, fs); - values[Traits::nonWettingPhaseIdx] = krn(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = krw(params, fs); + values[Traits::nonWettingPhaseIdx] = krn(params, fs); } /*! * \brief The capillary pressure-saturation curve */ - template - static Scalar pcnw(const Params ¶ms, const FluidState &fs) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) { - Scalar Sw = fs.saturation(Traits::wettingPhaseIdx); + typedef MathToolbox FsToolbox; + const auto& Sw = + FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); return twoPhaseSatPcnw(params, Sw); } @@ -130,117 +137,94 @@ public: /*! * \brief The saturation-capillary pressure curve */ - static Scalar twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) { return eval_(params.SwSamples(), params.pcnwSamples(), Sw); } /*! * \brief The saturation-capillary pressure curve */ - template - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) { OPM_THROW(std::logic_error, "Not implemented: Sw()"); } - static Scalar twoPhaseSatSw(const Params ¶ms, Scalar pC) + template + static Evaluation twoPhaseSatSw(const Params ¶ms, const Evaluation& pC) { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); } /*! * \brief Calculate the non-wetting phase saturations depending on * the phase pressures. */ - template - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return 1 - Sw(params, fs); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) + { return 1 - Sw(params, fs); } - static Scalar twoPhaseSatSn(const Params ¶ms, Scalar pC) + template + static Evaluation twoPhaseSatSn(const Params ¶ms, const Evaluation& pC) { return 1 - twoPhaseSatSw(params, pC); } - /*! - * \brief The partial derivative of the capillary pressure with - * regard to the saturation - */ - template - static Scalar dPcnw_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDPcnw_dSw(const Params ¶ms, 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 - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatKrw(const Params ¶ms, 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 - static Scalar dKrw_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDKrw_dSw(const Params ¶ms, Scalar Sw) + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) { - return evalDeriv_(params.SwSamples(), - params.krwSamples(), - Sw); + typedef MathToolbox FsToolbox; + const auto& Sw = + FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + + return twoPhaseSatKrw(params, Sw); + } + + template + static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation Sw) + { + typedef MathToolbox 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 - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); } - - static Scalar twoPhaseSatKrn(const Params ¶ms, Scalar Sw) + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) { - return std::max(0.0, std::min(1.0, eval_(params.SwSamples(), - params.krnSamples(), - Sw))); + typedef MathToolbox FsToolbox; + const auto& Sw = + FsToolbox::template toLhs(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 - static Scalar dKrn_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDKrn_dSw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) { - if (Sw < params.SwSamples().front() || Sw > params.SwSamples().back()) - return 0.0; - return evalDeriv_(params.SwSamples(), - params.krnSamples(), - Sw); + typedef MathToolbox 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 + 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 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 + static Evaluation evalDeriv_(const ValueVector &xValues, + const ValueVector &yValues, + const Evaluation& x) { - int segIdx = findSegmentIndex_(xValues, x); + typedef MathToolbox 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) diff --git a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp index 761ffdc52..2c043875a 100644 --- a/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp +++ b/opm/material/fluidmatrixinteractions/PiecewiseLinearTwoPhaseMaterialParams.hpp @@ -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 + 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 + 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 + 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 diff --git a/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp b/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp index 71616c5ae..26796bc7a 100644 --- a/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp +++ b/opm/material/fluidmatrixinteractions/RegularizedBrooksCorey.hpp @@ -115,8 +115,10 @@ public: template static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) { + typedef typename std::remove_reference::type Evaluation; + values[Traits::wettingPhaseIdx] = 0.0; // reference phase - values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); + values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); } /*! @@ -126,7 +128,9 @@ public: template static void saturations(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = Sw(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = Sw(params, fs); values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx]; } @@ -143,8 +147,10 @@ public: template static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = krw(params, fs); - values[Traits::nonWettingPhaseIdx] = krn(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = krw(params, fs); + values[Traits::nonWettingPhaseIdx] = krn(params, fs); } /*! @@ -171,13 +177,19 @@ public: * * \sa BrooksCorey::pcnw */ - template - static Scalar pcnw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) { - const Scalar Sthres = params.thresholdSw(); + typedef MathToolbox FsToolbox; + + const auto& Sw = FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + return twoPhaseSatPcnw(params, Sw); + } + + template + static Evaluation twoPhaseSatPcnw(const Params ¶ms, 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 - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) { - Scalar pcnw = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx); - return twoPhaseSatSw(params, pcnw); + typedef MathToolbox FsToolbox; + + const Evaluation& pC = + FsToolbox::template toLhs(fs.pressure(Traits::nonWettingPhaseIdx)) + - FsToolbox::template toLhs(fs.pressure(Traits::wettingPhaseIdx)); + return twoPhaseSatSw(params, pC); } - static Scalar twoPhaseSatSw(const Params ¶ms, Scalar pcnw) + template + static Evaluation twoPhaseSatSw(const Params ¶ms, 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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return 1 - Sw(params, fs); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) + { return 1 - Sw(params, fs); } - static Scalar twoPhaseSatSn(const Params ¶ms, Scalar pcnw) + template + static Evaluation twoPhaseSatSn(const Params ¶ms, 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 ¶ms, 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 ¶ms, 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 - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatKrw(const Params ¶ms, Scalar Sw) + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) { - if (Sw <= 0.0) - return 0.0; - else if (Sw >= 1.0) - return 1.0; + typedef MathToolbox FsToolbox; - return BrooksCorey::twoPhaseSatKrw(params, Sw); + const auto& Sw = FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + return twoPhaseSatKrw(params, Sw); } - static Scalar twoPhaseSatDKrw_dSw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) { - if (Sw <= 0.0 || Sw >= 1.0) - return 0.0; + typedef MathToolbox 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 - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); } - - static Scalar twoPhaseSatKrn(const Params ¶ms, Scalar Sw) + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) { + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = + 1.0 - FsToolbox::template toLhs(fs.saturation(Traits::nonWettingPhaseIdx)); + return twoPhaseSatKrn(params, Sw); + } + + template + static Evaluation twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox 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 ¶ms, Scalar Sw) - { - if (Sw <= 0.0 || Sw >= 1.0) - return 0.0; - - return BrooksCorey::twoPhaseSatDKrn_dSw(params, Sw); - } - }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/RegularizedBrooksCoreyParams.hpp b/opm/material/fluidmatrixinteractions/RegularizedBrooksCoreyParams.hpp index 86a90679e..5e84f42b0 100644 --- a/opm/material/fluidmatrixinteractions/RegularizedBrooksCoreyParams.hpp +++ b/opm/material/fluidmatrixinteractions/RegularizedBrooksCoreyParams.hpp @@ -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 + #include namespace Opm { @@ -38,6 +41,7 @@ template class RegularizedBrooksCoreyParams : public Opm::BrooksCoreyParams { typedef Opm::BrooksCoreyParams BrooksCoreyParams; + typedef Opm::BrooksCorey 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 diff --git a/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp b/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp index 83dcb09dd..b7d5f4dd6 100644 --- a/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp +++ b/opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp @@ -113,8 +113,10 @@ public: template static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) { + typedef typename std::remove_reference::type Evaluation; + values[Traits::wettingPhaseIdx] = 0.0; // reference phase - values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); + values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); } /*! @@ -124,7 +126,9 @@ public: template static void saturations(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = Sw(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = Sw(params, fs); values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx]; } @@ -135,8 +139,10 @@ public: template static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = krw(params, fs); - values[Traits::nonWettingPhaseIdx] = krn(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = krw(params, fs); + values[Traits::nonWettingPhaseIdx] = krn(params, fs); } /*! @@ -151,11 +157,17 @@ public: * * \copydetails VanGenuchten::pC() */ - template - static Scalar pcnw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatPcnw(params, fs.saturation(Traits::wettingPhaseIdx)); } + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) + { + typedef MathToolbox FsToolbox; - static Scalar twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + const auto& Sw = FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + return twoPhaseSatPcnw(params, Sw); + } + + template + static Evaluation twoPhaseSatPcnw(const Params ¶ms, 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 - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) { - Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx); + typedef MathToolbox FsToolbox; + + const Evaluation& pC = + FsToolbox::template toLhs(fs.pressure(Traits::nonWettingPhaseIdx)) + - FsToolbox::template toLhs(fs.pressure(Traits::wettingPhaseIdx)); return twoPhaseSatSw(params, pC); } - static Scalar twoPhaseSatSw(const Params ¶ms, Scalar pC) + template + static Evaluation twoPhaseSatSw(const Params ¶ms, 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& 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(/*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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return 1 - Sw(params, fs); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) + { return 1 - Sw(params, fs); } - static Scalar twoPhaseSatSn(const Params ¶ms, Scalar pc) + template + static Evaluation twoPhaseSatSn(const Params ¶ms, 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 - static Scalar dPcnw_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDPcnw_dSw(const Params ¶ms, 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 - static Scalar dSw_dpC(const Params ¶ms, 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 - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatKrw(const Params ¶ms, Scalar Sw) + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) { - // regularize - if (Sw < 0) - return 0; - else if (Sw > 1) - return 1; + typedef MathToolbox FsToolbox; - return VanGenuchten::twoPhaseSatKrw(params, Sw); + const auto& Sw = FsToolbox::template toLhs(fs.saturation(Traits::wettingPhaseIdx)); + return twoPhaseSatKrw(params, Sw); } - static Scalar twoPhaseSatDKrw_dSw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) { - if (Sw <= 0.0 || Sw >= 1.0) - return 0.0; + typedef MathToolbox 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 - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); } - - static Scalar twoPhaseSatKrn(const Params ¶ms, Scalar Sw) + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) { - // regularize - if (Sw <= 0) - return 1; - else if (Sw >= 1) - return 0; + typedef MathToolbox FsToolbox; - return VanGenuchten::twoPhaseSatKrn(params, Sw); + const Evaluation& Sw = + 1.0 - FsToolbox::template toLhs(fs.saturation(Traits::nonWettingPhaseIdx)); + return twoPhaseSatKrn(params, Sw); } - static Scalar twoPhaseSatDKrn_dSw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) { - if (Sw <= 0.0 || Sw >= 1.0) - return 0.0; + typedef MathToolbox 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); } }; diff --git a/opm/material/fluidmatrixinteractions/RegularizedVanGenuchtenParams.hpp b/opm/material/fluidmatrixinteractions/RegularizedVanGenuchtenParams.hpp index c30e1fa41..b1e09ee4f 100644 --- a/opm/material/fluidmatrixinteractions/RegularizedVanGenuchtenParams.hpp +++ b/opm/material/fluidmatrixinteractions/RegularizedVanGenuchtenParams.hpp @@ -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_; diff --git a/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp b/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp index 273d8cc45..b26b150c8 100644 --- a/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp +++ b/opm/material/fluidmatrixinteractions/SplineTwoPhaseMaterial.hpp @@ -88,10 +88,12 @@ public: * \brief The capillary pressure-saturation curve. */ template - static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) + static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fluidState) { + typedef typename std::remove_reference::type Evaluation; + values[Traits::wettingPhaseIdx] = 0.0; // reference phase - values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); + values[Traits::nonWettingPhaseIdx] = pcnw(params, fluidState); } /*! @@ -99,26 +101,31 @@ public: * pressure differences. */ template - static void saturations(Container &values, const Params ¶ms, const FluidState &fs) + static void saturations(Container &values, const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not implemented: saturations()"); } /*! * \brief The relative permeabilities */ template - static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) + static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fluidState) { - values[Traits::wettingPhaseIdx] = krw(params, fs); - values[Traits::nonWettingPhaseIdx] = krn(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = krw(params, fluidState); + values[Traits::nonWettingPhaseIdx] = krn(params, fluidState); } /*! * \brief The capillary pressure-saturation curve */ - template - static Scalar pcnw(const Params ¶ms, const FluidState &fs) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fluidState) { - Scalar Sw = fs.saturation(Traits::wettingPhaseIdx); + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = + FsToolbox::template toLhs(fluidState.saturation(Traits::wettingPhaseIdx)); return twoPhaseSatPcnw(params, Sw); } @@ -126,89 +133,78 @@ public: /*! * \brief The saturation-capillary pressure curve */ - static Scalar twoPhaseSatPcnw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) { return params.pcnwSpline().eval(Sw, /*extrapolate=*/true); } /*! * \brief The saturation-capillary pressure curve */ - template - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not implemented: Sw()"); } - static Scalar twoPhaseSatSw(const Params ¶ms, Scalar pC) + template + static Evaluation twoPhaseSatSw(const Params ¶ms, const Evaluation& pC) { OPM_THROW(std::logic_error, "Not implemented: twoPhaseSatSw()"); } /*! * \brief Calculate the non-wetting phase saturations depending on * the phase pressures. */ - template - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return 1 - Sw(params, fs); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fluidState) + { return 1 - Sw(params, fluidState); } - static Scalar twoPhaseSatSn(const Params ¶ms, Scalar pC) + template + static Evaluation twoPhaseSatSn(const Params ¶ms, const Evaluation& pC) { return 1 - twoPhaseSatSw(params, pC); } - /*! - * \brief The partial derivative of the capillary pressure with - * regard to the saturation - */ - template - static Scalar dPcnw_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDPcnw_dSw(const Params ¶ms, 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 - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); } + template + static Evaluation krw(const Params ¶ms, const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; - static Scalar twoPhaseSatKrw(const Params ¶ms, Scalar Sw) - { return std::max(0.0, std::min(1.0, params.krwSpline().eval(Sw, /*extrapolate=*/true))); } + const Evaluation& Sw = + FsToolbox::template toLhs(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 - static Scalar dKrw_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } + return twoPhaseSatKrw(params, Sw); + } - static Scalar twoPhaseSatDKrw_dSw(const Params ¶ms, Scalar Sw) - { return params.krwSpline().evalDerivative(Sw, /*extrapolate=*/true); } + template + static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox 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 - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); } + template + static Evaluation krn(const Params ¶ms, const FluidState &fluidState) + { + typedef MathToolbox FsToolbox; - static Scalar twoPhaseSatKrn(const Params ¶ms, Scalar Sw) - { return std::max(0.0, std::min(1.0, params.krnSpline().eval(Sw, /*extrapolate=*/true))); } + const Evaluation& Sn = + FsToolbox::template toLhs(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 - static Scalar dKrn_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } + return twoPhaseSatKrn(params, 1.0 - Sn); + } - static Scalar twoPhaseSatDKrn_dSw(const Params ¶ms, Scalar Sw) - { return params.krnSpline().evalDerivative(Sw, /*extrapolate=*/true); } + template + static Evaluation twoPhaseSatKrn(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox Toolbox; + + return Toolbox::max(0.0, Toolbox::min(1.0, params.krnSpline().eval(Sw, /*extrapolate=*/true))); + } }; } // namespace Opm diff --git a/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp b/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp index 9ebcd343f..f0cce424a 100644 --- a/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp +++ b/opm/material/fluidmatrixinteractions/ThreePhaseParkerVanGenuchten.hpp @@ -27,6 +27,8 @@ #include "ThreePhaseParkerVanGenuchtenParams.hpp" +#include + #include namespace Opm { @@ -99,9 +101,11 @@ public: const Params ¶ms, const FluidState &fluidState) { - values[gasPhaseIdx] = pcgn(params, fluidState); + typedef typename std::remove_reference::type Evaluation; + + values[gasPhaseIdx] = pcgn(params, fluidState); values[nonWettingPhaseIdx] = 0; - values[wettingPhaseIdx] = - pcnw(params, fluidState); + values[wettingPhaseIdx] = - pcnw(params, fluidState); } /*! @@ -113,18 +117,20 @@ public: * p_{c,gn} = p_g - p_n * \f] */ - template - static Scalar pcgn(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation pcgn(const Params ¶ms, const FluidState &fluidState) { + typedef MathToolbox FsToolbox; + typedef MathToolbox 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(fluidState.saturation(wettingPhaseIdx)) + + FsToolbox::template toLhs(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 - static Scalar pcnw(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fluidState) { - Scalar Sw = fluidState.saturation(wettingPhaseIdx); - Scalar Se = (Sw-params.Swr())/(1.-params.Snr()); + typedef MathToolbox FsToolbox; + typedef MathToolbox Toolbox; + + const Evaluation& Sw = + FsToolbox::template toLhs(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 - static Scalar Sg(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation Sg(const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not implemented: Sg()"); } /*! * \brief The saturation of the non-wetting (i.e., oil) phase. */ - template - static Scalar Sn(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation Sn(const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not implemented: Sn()"); } /*! * \brief The saturation of the wetting (i.e., water) phase. */ - template - static Scalar Sw(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fluidState) { OPM_THROW(std::logic_error, "Not implemented: Sw()"); } /*! @@ -251,9 +257,11 @@ public: const Params ¶ms, const FluidState &fluidState) { - values[wettingPhaseIdx] = krw(params, fluidState); - values[nonWettingPhaseIdx] = krn(params, fluidState); - values[gasPhaseIdx] = krg(params, fluidState); + typedef typename std::remove_reference::type Evaluation; + + values[wettingPhaseIdx] = krw(params, fluidState); + values[nonWettingPhaseIdx] = krn(params, fluidState); + values[gasPhaseIdx] = krg(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 - static Scalar krw(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation krw(const Params ¶ms, const FluidState &fluidState) { + typedef MathToolbox FsToolbox; + typedef MathToolbox Toolbox; + const Evaluation& Sw = + FsToolbox::template toLhs(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 - static Scalar krn(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation krn(const Params ¶ms, 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 FsToolbox; + typedef MathToolbox Toolbox; + + const Evaluation& Sn = + FsToolbox::template toLhs(fluidState.saturation(nonWettingPhaseIdx)); + const Evaluation& Sw = + FsToolbox::template toLhs(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 - static Scalar krg(const Params ¶ms, - const FluidState &fluidState) + template + static Evaluation krg(const Params ¶ms, const FluidState &fluidState) { - Scalar Sg = fluidState.saturation(gasPhaseIdx); - Scalar Se = std::min(((1-Sg) - params.Sgr()) / (1 - params.Sgr()), 1.); + typedef MathToolbox FsToolbox; + typedef MathToolbox Toolbox; + + const Evaluation& Sg = + FsToolbox::template toLhs(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 diff --git a/opm/material/fluidmatrixinteractions/VanGenuchten.hpp b/opm/material/fluidmatrixinteractions/VanGenuchten.hpp index 3e63c9f13..161fee5ab 100644 --- a/opm/material/fluidmatrixinteractions/VanGenuchten.hpp +++ b/opm/material/fluidmatrixinteractions/VanGenuchten.hpp @@ -26,6 +26,8 @@ #include "VanGenuchtenParams.hpp" +#include + #include #include #include @@ -108,8 +110,10 @@ public: template static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) { + typedef typename std::remove_reference::type Evaluation; + values[Traits::wettingPhaseIdx] = 0.0; // reference phase - values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); + values[Traits::nonWettingPhaseIdx] = pcnw(params, fs); } /*! @@ -119,7 +123,9 @@ public: template static void saturations(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = Sw(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = Sw(params, fs); values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx]; } @@ -136,8 +142,10 @@ public: template static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) { - values[Traits::wettingPhaseIdx] = krw(params, fs); - values[Traits::nonWettingPhaseIdx] = krn(params, fs); + typedef typename std::remove_reference::type Evaluation; + + values[Traits::wettingPhaseIdx] = krw(params, fs); + values[Traits::nonWettingPhaseIdx] = krn(params, fs); } /*! @@ -154,10 +162,14 @@ public: * \param fs The fluid state for which the capillary pressure * ought to be calculated */ - template - static Scalar pcnw(const Params ¶ms, const FluidState &fs) + template + static Evaluation pcnw(const Params ¶ms, const FluidState &fs) { - Scalar Sw = fs.saturation(Traits::wettingPhaseIdx); + typedef MathToolbox FsToolbox; + + const Evaluation& Sw = + FsToolbox::template toLhs(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 ¶ms, Scalar Sw) - { return std::pow(std::pow(Sw, -1.0/params.vgM()) - 1, 1.0/params.vgN())/params.vgAlpha(); } + template + static Evaluation twoPhaseSatPcnw(const Params ¶ms, const Evaluation& Sw) + { + typedef MathToolbox 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 - static Scalar Sw(const Params ¶ms, const FluidState &fs) + template + static Evaluation Sw(const Params ¶ms, const FluidState &fs) { - Scalar pC = fs.pressure(Traits::nonWettingPhaseIdx) - fs.pressure(Traits::wettingPhaseIdx); + typedef MathToolbox FsToolbox; + + Evaluation pC = + FsToolbox::template toLhs(fs.pressure(Traits::nonWettingPhaseIdx)) + - FsToolbox::template toLhs(fs.pressure(Traits::wettingPhaseIdx)); return twoPhaseSatSw(params, pC); } - static Scalar twoPhaseSatSw(const Params ¶ms, Scalar pC) + template + static Evaluation twoPhaseSatSw(const Params ¶ms, const Evaluation& pC) { + typedef MathToolbox 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 - static Scalar Sn(const Params ¶ms, const FluidState &fs) - { return 1 - Sw(params, fs); } + template + static Evaluation Sn(const Params ¶ms, const FluidState &fs) + { return 1 - Sw(params, fs); } - static Scalar twoPhaseSatSn(const Params ¶ms, Scalar pC) + template + static Evaluation twoPhaseSatSn(const Params ¶ms, 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 - static Scalar dPcnw_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDPcnw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDPcnw_dSw(const Params ¶ms, 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 - static Scalar krw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatKrw(const Params ¶ms, Scalar Sw) + template + static Evaluation krw(const Params ¶ms, const FluidState &fs) { - assert(0 <= Sw && Sw <= 1); + typedef MathToolbox 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(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 - static Scalar dKrw_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDkrw_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDKrw_dSw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatKrw(const Params ¶ms, const Evaluation& Sw) { - assert(0 <= Sw && Sw <= 1); + typedef MathToolbox 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 - static Scalar krn(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatKrn(params, 1.0 - fs.saturation(Traits::nonWettingPhaseIdx)); } - - static Scalar twoPhaseSatKrn(const Params ¶ms, Scalar Sw) + template + static Evaluation krn(const Params ¶ms, const FluidState &fs) { - assert(0 <= Sw && Sw <= 1); + typedef MathToolbox 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(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 - static Scalar dKrn_dSw(const Params ¶ms, const FluidState &fs) - { return twoPhaseSatDkrn_dSw(params, fs.saturation(Traits::wettingPhaseIdx)); } - - static Scalar twoPhaseSatDKrn_dSw(const Params ¶ms, Scalar Sw) + template + static Evaluation twoPhaseSatKrn(const Params ¶ms, Evaluation Sw) { + typedef MathToolbox 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 diff --git a/tests/test_fluidmatrixinteractions.cpp b/tests/test_fluidmatrixinteractions.cpp index 0132b7468..d0715c603 100644 --- a/tests/test_fluidmatrixinteractions.cpp +++ b/tests/test_fluidmatrixinteractions.cpp @@ -25,6 +25,9 @@ */ #include "config.h" +// include the local AD framwork +#include + // include all capillary pressure laws #include #include @@ -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(params, fs); + v = MaterialLaw::template Sw(params, fs); + v = MaterialLaw::template Sn(params, fs); + v = MaterialLaw::template krw(params, fs); + v = MaterialLaw::template krn(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(params, fs); + v = MaterialLaw::template Sw(params, fs); + v = MaterialLaw::template Sn(params, fs); + v = MaterialLaw::template Sg(params, fs); + v = MaterialLaw::template krw(params, fs); + v = MaterialLaw::template krn(params, fs); + v = MaterialLaw::template krg(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 TwoPhaseFluidState; - typedef Opm::ImmiscibleFluidState ThreePhaseFluidState; + typedef Opm::LocalAd::Evaluation Evaluation; + typedef Opm::ImmiscibleFluidState TwoPhaseFluidState; + typedef Opm::ImmiscibleFluidState ThreePhaseFluidState; MyMpiHelper mpiHelper(argc, argv);