From 26c85539027664874af13827f450969262d4eac0 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 21 May 2015 15:33:16 +0200 Subject: [PATCH] make all fluid systems local-AD aware also, adapt the unit test make sure they can be used synchronously with function evaluations as well as scalars. --- opm/material/fluidsystems/BaseFluidSystem.hpp | 78 ++--- .../fluidsystems/BlackOilFluidSystem.hpp | 239 +++++++++------ .../fluidsystems/BrineCO2FluidSystem.hpp | 274 +++++++++-------- opm/material/fluidsystems/GasPhase.hpp | 40 ++- .../fluidsystems/H2OAirFluidSystem.hpp | 158 +++++----- .../H2OAirMesityleneFluidSystem.hpp | 204 ++++++------- .../fluidsystems/H2OAirXyleneFluidSystem.hpp | 203 ++++++------- .../fluidsystems/H2ON2FluidSystem.hpp | 276 +++++++++--------- .../H2ON2LiquidPhaseFluidSystem.hpp | 124 ++++---- opm/material/fluidsystems/LiquidPhase.hpp | 24 +- .../fluidsystems/SinglePhaseFluidSystem.hpp | 98 ++++--- opm/material/fluidsystems/Spe5FluidSystem.hpp | 12 +- .../TwoPhaseImmiscibleFluidSystem.hpp | 88 +++--- .../ConstantCompressibilityOilPvt.hpp | 111 ++++--- .../ConstantCompressibilityWaterPvt.hpp | 49 ++-- .../fluidsystems/blackoilpvt/DeadOilPvt.hpp | 113 ++++--- .../fluidsystems/blackoilpvt/DryGasPvt.hpp | 106 ++++--- .../blackoilpvt/GasPvtInterface.hpp | 207 ++++++++++++- .../fluidsystems/blackoilpvt/LiveOilPvt.hpp | 201 +++++++------ .../blackoilpvt/OilPvtInterface.hpp | 251 +++++++++++++++- .../blackoilpvt/WaterPvtInterface.hpp | 114 +++++++- .../fluidsystems/blackoilpvt/WetGasPvt.hpp | 158 +++++----- tests/checkFluidSystem.hpp | 39 ++- tests/test_fluidsystems.cpp | 135 +++++---- 24 files changed, 2100 insertions(+), 1202 deletions(-) diff --git a/opm/material/fluidsystems/BaseFluidSystem.hpp b/opm/material/fluidsystems/BaseFluidSystem.hpp index 826627bea..3c6f382b0 100644 --- a/opm/material/fluidsystems/BaseFluidSystem.hpp +++ b/opm/material/fluidsystems/BaseFluidSystem.hpp @@ -65,7 +65,7 @@ public: static char *phaseName(int phaseIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a phaseName() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a phaseName() method!"); } /*! @@ -76,7 +76,7 @@ public: static bool isLiquid(int phaseIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isLiquid() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isLiquid() method!"); } /*! @@ -96,7 +96,7 @@ public: static bool isIdealMixture(int phaseIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isIdealMixture() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isIdealMixture() method!"); } /*! @@ -111,7 +111,7 @@ public: static bool isCompressible(int phaseIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isCompressible() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isCompressible() method!"); } /*! @@ -123,7 +123,7 @@ public: static bool isIdealGas(int phaseIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isIdealGas() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a isIdealGas() method!"); } /*! @@ -134,7 +134,7 @@ public: static const char *componentName(int compIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a componentName() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a componentName() method!"); } /*! @@ -145,7 +145,7 @@ public: static Scalar molarMass(int compIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a molarMass() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a molarMass() method!"); } /*! @@ -160,13 +160,13 @@ public: * \copydoc Doxygen::fluidSystemBaseParams * \copydoc Doxygen::phaseIdxParam */ - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { OPM_THROW(std::runtime_error, - "Not implemented: The fluid system '" << Opm::className() << "' does not provide a density() method!"); + "Not implemented: The fluid system '" << Opm::className() << "' does not provide a density() method!"); } /*! @@ -183,11 +183,11 @@ public: * \copydoc Doxygen::phaseIdxParam * \copydoc Doxygen::compIdxParam */ - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a fugacityCoefficient() method!"); } @@ -198,10 +198,10 @@ public: * \copydoc Doxygen::fluidSystemBaseParams * \copydoc Doxygen::phaseIdxParam */ - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a viscosity() method!"); } @@ -223,13 +223,13 @@ public: * \copydoc Doxygen::phaseIdxParam * \copydoc Doxygen::compIdxParam */ - template - static Scalar diffusionCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval diffusionCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { - OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a diffusionCoefficient() method!"); + OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a diffusionCoefficient() method!"); } /*! @@ -239,10 +239,10 @@ public: * \copydoc Doxygen::fluidSystemBaseParams * \copydoc Doxygen::phaseIdxParam */ - template - static Scalar enthalpy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide an enthalpy() method!"); } @@ -253,10 +253,10 @@ public: * \copydoc Doxygen::fluidSystemBaseParams * \copydoc Doxygen::phaseIdxParam */ - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a thermalConductivity() method!"); } @@ -267,10 +267,10 @@ public: * \copydoc Doxygen::fluidSystemBaseParams * \copydoc Doxygen::phaseIdxParam */ - template - static Scalar heatCapacity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval heatCapacity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className() << "' does not provide a heatCapacity() method!"); } diff --git a/opm/material/fluidsystems/BlackOilFluidSystem.hpp b/opm/material/fluidsystems/BlackOilFluidSystem.hpp index 94359c35f..e1de044a4 100644 --- a/opm/material/fluidsystems/BlackOilFluidSystem.hpp +++ b/opm/material/fluidsystems/BlackOilFluidSystem.hpp @@ -38,27 +38,17 @@ #include namespace Opm { -template -class OilPvtInterface; - -template -class GasPvtInterface; - -template -class WaterPvtInterface; - namespace FluidSystems { - /*! * \brief A fluid system which uses the black-oil parameters * to calculate termodynamically meaningful quantities. */ -template -class BlackOil : public BaseFluidSystem > +template +class BlackOil : public BaseFluidSystem > { - typedef Opm::GasPvtInterface GasPvtInterface; - typedef Opm::OilPvtInterface OilPvtInterface; - typedef Opm::WaterPvtInterface WaterPvtInterface; + typedef Opm::GasPvtInterface GasPvtInterface; + typedef Opm::OilPvtInterface OilPvtInterface; + typedef Opm::WaterPvtInterface WaterPvtInterface; public: //! \copydoc BaseFluidSystem::ParameterCache @@ -286,26 +276,29 @@ public: * thermodynamic relations ****************************************/ //! \copydoc BaseFluidSystem::density - template - static Scalar density(const FluidState &fluidState, - ParameterCache ¶mCache, - const int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + ParameterCache ¶mCache, + const int phaseIdx) { assert(0 <= phaseIdx && phaseIdx <= numPhases); - Scalar p = fluidState.pressure(phaseIdx); - Scalar T = fluidState.temperature(phaseIdx); + typedef typename FluidState::Scalar FsEval; + typedef Opm::MathToolbox FsToolbox; + + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); int regionIdx = paramCache.regionIndex(); switch (phaseIdx) { - case waterPhaseIdx: return waterDensity(T, p, regionIdx); + case waterPhaseIdx: return waterDensity(T, p, regionIdx); case gasPhaseIdx: { - Scalar XgO = fluidState.massFraction(gasPhaseIdx, oilCompIdx); - return gasDensity(T, p, XgO, regionIdx); + const auto& XgO = FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, oilCompIdx)); + return gasDensity(T, p, XgO, regionIdx); } case oilPhaseIdx: { - Scalar XoG = fluidState.massFraction(oilPhaseIdx, gasCompIdx); - return oilDensity(T, p, XoG, regionIdx); + const auto& XoG = FsToolbox::template toLhs(fluidState.massFraction(oilPhaseIdx, gasCompIdx)); + return oilDensity(T, p, XoG, regionIdx); } } @@ -313,49 +306,53 @@ public: } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { assert(0 <= phaseIdx && phaseIdx <= numPhases); assert(0 <= compIdx && compIdx <= numComponents); - Scalar p = fluidState.pressure(phaseIdx); - Scalar T = fluidState.temperature(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); int regionIdx = paramCache.regionIndex(); switch (phaseIdx) { - case waterPhaseIdx: return fugCoefficientInWater(compIdx, T, p, regionIdx); - case gasPhaseIdx: return fugCoefficientInGas(compIdx, T, p, regionIdx); - case oilPhaseIdx: return fugCoefficientInOil(compIdx, T, p, regionIdx); + case waterPhaseIdx: return fugCoefficientInWater(compIdx, T, p, regionIdx); + case gasPhaseIdx: return fugCoefficientInGas(compIdx, T, p, regionIdx); + case oilPhaseIdx: return fugCoefficientInOil(compIdx, T, p, regionIdx); } OPM_THROW(std::logic_error, "Unhandled phase or component index"); } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { assert(0 <= phaseIdx && phaseIdx <= numPhases); - Scalar p = fluidState.pressure(phaseIdx); - Scalar T = fluidState.temperature(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); int regionIdx = paramCache.regionIndex(); switch (phaseIdx) { case oilPhaseIdx: { - Scalar XoG = fluidState.massFraction(oilPhaseIdx, gasCompIdx); + const auto& XoG = FsToolbox::template toLhs(fluidState.massFraction(oilPhaseIdx, gasCompIdx)); return oilPvt_->viscosity(regionIdx, T, p, XoG); } case waterPhaseIdx: return waterPvt_->viscosity(regionIdx, T, p); case gasPhaseIdx: { - Scalar XgO = fluidState.massFraction(gasPhaseIdx, oilCompIdx); + const auto& XgO = FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, oilCompIdx)); return gasPvt_->viscosity(regionIdx, T, p, XgO); } } @@ -394,14 +391,15 @@ public: * * \param pressure The pressure of interest [Pa] */ - static Scalar saturatedOilFormationVolumeFactor(Scalar temperature, - Scalar pressure, - int regionIdx) + template + static LhsEval saturatedOilFormationVolumeFactor(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { Valgrind::CheckDefined(pressure); // calculate the mass fractions of gas and oil - Scalar XoG = saturatedOilGasMassFraction(temperature, pressure, regionIdx); + const auto& XoG = saturatedOilGasMassFraction(temperature, pressure, regionIdx); // ATTENTION: XoG is represented by the _first_ axis! return oilFormationVolumeFactor(temperature, pressure, XoG, regionIdx); @@ -410,7 +408,10 @@ public: /*! * \brief Return the formation volume factor of water. */ - static Scalar waterFormationVolumeFactor(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval waterFormationVolumeFactor(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return waterPvt_->formationVolumeFactor(regionIdx, temperature, pressure); } /*! @@ -418,7 +419,10 @@ public: * * \param pressure The pressure of interest [Pa] */ - static Scalar gasDissolutionFactor(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval gasDissolutionFactor(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return oilPvt_->gasDissolutionFactor(regionIdx, temperature, pressure); } /*! @@ -426,7 +430,10 @@ public: * * \param pressure The pressure of interest [Pa] */ - static Scalar oilVaporizationFactor(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval oilVaporizationFactor(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return gasPvt_->oilVaporizationFactor(regionIdx, temperature, pressure); } /*! @@ -435,7 +442,11 @@ public: * \param compIdx The index of the component of interest * \param pressure The pressure of interest [Pa] */ - static Scalar fugCoefficientInWater(int compIdx, Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval fugCoefficientInWater(int compIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return waterPvt_->fugacityCoefficient(regionIdx, temperature, pressure, compIdx); } /*! @@ -444,7 +455,11 @@ public: * \param compIdx The index of the component of interest * \param pressure The pressure of interest [Pa] */ - static Scalar fugCoefficientInGas(int compIdx, Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval fugCoefficientInGas(int compIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return gasPvt_->fugacityCoefficient(regionIdx, temperature, pressure, compIdx); } /*! @@ -453,7 +468,11 @@ public: * \param compIdx The index of the component of interest * \param pressure The pressure of interest [Pa] */ - static Scalar fugCoefficientInOil(int compIdx, Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval fugCoefficientInOil(int compIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return oilPvt_->fugacityCoefficient(regionIdx, temperature, pressure, compIdx); } /*! @@ -462,75 +481,109 @@ public: * * \param XoG The mass fraction of the gas component in the oil phase [-] */ - static Scalar oilSaturationPressure(Scalar temperature, Scalar XoG, int regionIdx) + template + static LhsEval oilSaturationPressure(const LhsEval& temperature, + const LhsEval& XoG, + int regionIdx) { return oilPvt_->oilSaturationPressure(regionIdx, temperature, XoG); } /*! * \brief The maximum mass fraction of the gas component in the oil phase. */ - static Scalar saturatedOilGasMassFraction(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval saturatedOilGasMassFraction(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return oilPvt_->saturatedOilGasMassFraction(regionIdx, temperature, pressure); } /*! * \brief The maximum mole fraction of the gas component in the oil phase. */ - static Scalar saturatedOilGasMoleFraction(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval saturatedOilGasMoleFraction(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return oilPvt_->saturatedOilGasMoleFraction(regionIdx, temperature, pressure); } /*! * \brief The maximum mass fraction of the oil component in the gas phase. */ - static Scalar saturatedGasOilMassFraction(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval saturatedGasOilMassFraction(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return gasPvt_->saturatedGasOilMassFraction(regionIdx, temperature, pressure); } /*! * \brief The maximum mole fraction of the oil component in the gas phase. */ - static Scalar saturatedGasOilMoleFraction(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval saturatedGasOilMoleFraction(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return gasPvt_->saturatedGasOilMoleFraction(regionIdx, temperature, pressure); } /*! * \brief Return the normalized formation volume factor of (potentially) * under-saturated oil. */ - static Scalar oilFormationVolumeFactor(Scalar temperature, - Scalar pressure, - Scalar XoG, - int regionIdx) + template + static LhsEval oilFormationVolumeFactor(const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG, + int regionIdx) { return oilPvt_->formationVolumeFactor(regionIdx, temperature, pressure, XoG); } /*! * \brief Return the density of (potentially) under-saturated oil. */ - static Scalar oilDensity(Scalar temperature, Scalar pressure, Scalar XoG, int regionIdx) + template + static LhsEval oilDensity(const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG, + int regionIdx) { return oilPvt_->density(regionIdx, temperature, pressure, XoG); } /*! * \brief Return the density of gas-saturated oil. */ - static Scalar saturatedOilDensity(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval saturatedOilDensity(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { // mass fraction of gas-saturated oil - Scalar XoG = saturatedOilGasMassFraction(temperature, pressure, regionIdx); + const LhsEval& XoG = saturatedOilGasMassFraction(temperature, pressure, regionIdx); return oilPvt_->density(regionIdx, temperature, pressure, XoG); } /*! * \brief Return the formation volume factor of gas. */ - static Scalar gasFormationVolumeFactor(Scalar temperature, Scalar pressure, Scalar XgO, int regionIdx) + template + static LhsEval gasFormationVolumeFactor(const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO, + int regionIdx) { return gasPvt_->formationVolumeFactor(regionIdx, temperature, pressure, XgO); } /*! * \brief Return the density of dry gas. */ - static Scalar gasDensity(Scalar temperature, Scalar pressure, Scalar XgO, int regionIdx) + template + static LhsEval gasDensity(const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO, + int regionIdx) { return gasPvt_->density(regionIdx, temperature, pressure, XgO); } /*! * \brief Return the density of water. */ - static Scalar waterDensity(Scalar temperature, Scalar pressure, int regionIdx) + template + static LhsEval waterDensity(const LhsEval& temperature, + const LhsEval& pressure, + int regionIdx) { return waterPvt_->density(regionIdx, temperature, pressure); } private: @@ -540,9 +593,9 @@ private: referenceDensity_.resize(numRegions); } - static std::shared_ptr > gasPvt_; - static std::shared_ptr > oilPvt_; - static std::shared_ptr > waterPvt_; + static std::shared_ptr > gasPvt_; + static std::shared_ptr > oilPvt_; + static std::shared_ptr > waterPvt_; static bool enableDissolvedGas_; static bool enableVaporizedOil_; @@ -554,39 +607,39 @@ private: static std::vector > molarMass_; }; -template +template const Scalar -BlackOil::surfaceTemperature = 273.15 + 15.56; // [K] +BlackOil::surfaceTemperature = 273.15 + 15.56; // [K] -template +template const Scalar -BlackOil::surfacePressure = 101325.0; // [Pa] +BlackOil::surfacePressure = 101325.0; // [Pa] -template -bool BlackOil::enableDissolvedGas_; +template +bool BlackOil::enableDissolvedGas_; -template -bool BlackOil::enableVaporizedOil_; +template +bool BlackOil::enableVaporizedOil_; -template -std::shared_ptr > -BlackOil::oilPvt_; +template +std::shared_ptr > +BlackOil::oilPvt_; -template -std::shared_ptr > -BlackOil::gasPvt_; +template +std::shared_ptr > +BlackOil::gasPvt_; -template -std::shared_ptr > -BlackOil::waterPvt_; +template +std::shared_ptr > +BlackOil::waterPvt_; -template +template std::vector > -BlackOil::referenceDensity_; +BlackOil::referenceDensity_; -template +template std::vector > -BlackOil::molarMass_; +BlackOil::molarMass_; }} // namespace Opm, FluidSystems #endif diff --git a/opm/material/fluidsystems/BrineCO2FluidSystem.hpp b/opm/material/fluidsystems/BrineCO2FluidSystem.hpp index d039dff5c..37d78fdef 100644 --- a/opm/material/fluidsystems/BrineCO2FluidSystem.hpp +++ b/opm/material/fluidsystems/BrineCO2FluidSystem.hpp @@ -226,30 +226,33 @@ public: /*! * \copydoc BaseFluidSystem::density */ - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + typedef MathToolbox LhsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const LhsEval& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == liquidPhaseIdx) { // use normalized composition for to calculate the density // (the relations don't seem to take non-normalized // compositions too well...) - Scalar xlBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(liquidPhaseIdx, BrineIdx))); - Scalar xlCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(liquidPhaseIdx, CO2Idx))); - Scalar sumx = xlBrine + xlCO2; + LhsEval xlBrine = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, BrineIdx)))); + LhsEval xlCO2 = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, CO2Idx)))); + LhsEval sumx = xlBrine + xlCO2; xlBrine /= sumx; xlCO2 /= sumx; - Scalar result = liquidDensity_(temperature, - pressure, - xlBrine, - xlCO2); + LhsEval result = liquidDensity_(temperature, + pressure, + xlBrine, + xlCO2); Valgrind::CheckDefined(result); return result; @@ -260,16 +263,16 @@ public: // use normalized composition for to calculate the density // (the relations don't seem to take non-normalized // compositions too well...) - Scalar xgBrine = std::min(1.0, std::max(0.0, fluidState.moleFraction(gasPhaseIdx, BrineIdx))); - Scalar xgCO2 = std::min(1.0, std::max(0.0, fluidState.moleFraction(gasPhaseIdx, CO2Idx))); - Scalar sumx = xgBrine + xgCO2; + LhsEval xgBrine = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, BrineIdx)))); + LhsEval xgCO2 = LhsToolbox::min(1.0, LhsToolbox::max(0.0, FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, CO2Idx)))); + LhsEval sumx = xgBrine + xgCO2; xgBrine /= sumx; xgCO2 /= sumx; - Scalar result = gasDensity_(temperature, - pressure, - xgBrine, - xgCO2); + LhsEval result = gasDensity_(temperature, + pressure, + xgBrine, + xgCO2); Valgrind::CheckDefined(result); return result; } @@ -277,26 +280,28 @@ public: /*! * \copydoc BaseFluidSystem::viscosity */ - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const LhsEval& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == liquidPhaseIdx) { // assume pure brine for the liquid phase. TODO: viscosity // of mixture - Scalar result = Brine::liquidViscosity(temperature, pressure); + LhsEval result = Brine::liquidViscosity(temperature, pressure); Valgrind::CheckDefined(result); return result; } assert(phaseIdx == gasPhaseIdx); - Scalar result = CO2::gasViscosity(temperature, pressure); + LhsEval result = CO2::gasViscosity(temperature, pressure); Valgrind::CheckDefined(result); return result; } @@ -304,12 +309,15 @@ public: /*! * \copydoc BaseFluidSystem::fugacityCoefficient */ - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { + typedef MathToolbox FsToolbox; + typedef MathToolbox LhsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); @@ -317,18 +325,18 @@ public: // use the fugacity coefficients of an ideal gas. the // actual value of the fugacity is not relevant, as long // as the relative fluid compositions are observed, - return 1.0; + return LhsToolbox::createConstant(1.0); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const LhsEval& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); assert(temperature > 0); assert(pressure > 0); // calulate the equilibrium composition for the given // temperature and pressure. TODO: calculateMoleFractions() // could use some cleanup. - Scalar xlH2O, xgH2O; - Scalar xlCO2, xgCO2; + LhsEval xlH2O, xgH2O; + LhsEval xlCO2, xgCO2; BinaryCoeffBrineCO2::calculateMoleFractions(temperature, pressure, Brine_IAPWS::salinity, @@ -337,8 +345,8 @@ public: xgH2O); // normalize the phase compositions - xlCO2 = std::max(0.0, std::min(1.0, xlCO2)); - xgH2O = std::max(0.0, std::min(1.0, xgH2O)); + xlCO2 = LhsToolbox::max(0.0, LhsToolbox::min(1.0, xlCO2)); + xgH2O = LhsToolbox::max(0.0, LhsToolbox::min(1.0, xgH2O)); xlH2O = 1.0 - xlCO2; xgCO2 = 1.0 - xgH2O; @@ -358,14 +366,16 @@ public: /*! * \copydoc BaseFluidSystem::diffusionCoefficient */ - template - static Scalar diffusionCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval diffusionCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + typedef MathToolbox FsToolbox; + + const LhsEval& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == liquidPhaseIdx) return BinaryCoeffBrineCO2::liquidDiffCoeff(temperature, pressure); @@ -376,30 +386,33 @@ public: /*! * \copydoc BaseFluidSystem::enthalpy */ - template - static Scalar enthalpy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + typedef MathToolbox LhsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const LhsEval& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == liquidPhaseIdx) { - Scalar XlCO2 = fluidState.massFraction(phaseIdx, CO2Idx); - Scalar result = liquidEnthalpyBrineCO2_(temperature, - pressure, - Brine_IAPWS::salinity, - XlCO2); + const LhsEval& XlCO2 = FsToolbox::template toLhs(fluidState.massFraction(phaseIdx, CO2Idx)); + const LhsEval& result = liquidEnthalpyBrineCO2_(temperature, + pressure, + Brine_IAPWS::salinity, + XlCO2); Valgrind::CheckDefined(result); return result; } else { - Scalar XCO2 = fluidState.massFraction(gasPhaseIdx, CO2Idx); - Scalar XBrine = fluidState.massFraction(gasPhaseIdx, BrineIdx); + const LhsEval& XCO2 = FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, CO2Idx)); + const LhsEval& XBrine = FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, BrineIdx)); - Scalar result = 0; + LhsEval result = LhsToolbox::createConstant(0); result += XBrine * Brine::gasEnthalpy(temperature, pressure); result += XCO2 * CO2::gasEnthalpy(temperature, pressure); Valgrind::CheckDefined(result); @@ -410,17 +423,19 @@ public: /*! * \copydoc BaseFluidSystem::thermalConductivity */ - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox LhsToolbox; + // TODO way too simple! if (phaseIdx == liquidPhaseIdx) - return 0.6; // conductivity of water[W / (m K ) ] + return LhsToolbox::createConstant(0.6); // conductivity of water[W / (m K ) ] // gas phase - return 0.025; // conductivity of air [W / (m K ) ] + return LhsToolbox::createConstant(0.025); // conductivity of air [W / (m K ) ] } /*! @@ -435,32 +450,37 @@ public: * \param phaseIdx The index of the fluid phase to consider * \tparam FluidState the fluid state class */ - template - static Scalar heatCapacity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval heatCapacity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + + assert(0 <= phaseIdx && phaseIdx < numPhases); + + const LhsEval& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + if(phaseIdx == liquidPhaseIdx) - return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx)); + return H2O::liquidHeatCapacity(temperature, pressure); else - return CO2::gasHeatCapacity(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx)); + return CO2::gasHeatCapacity(temperature, pressure); } private: - static Scalar gasDensity_(Scalar T, - Scalar pg, - Scalar xgH2O, - Scalar xgCO2) + template + static LhsEval gasDensity_(const LhsEval& T, + const LhsEval& pg, + const LhsEval& xgH2O, + const LhsEval& xgCO2) { Valgrind::CheckDefined(T); Valgrind::CheckDefined(pg); Valgrind::CheckDefined(xgH2O); Valgrind::CheckDefined(xgCO2); - Scalar gasDensity = CO2::gasDensity(T, pg); - return gasDensity; + return CO2::gasDensity(T, pg); } /***********************************************************************/ @@ -469,10 +489,11 @@ private: /* rho_{b,CO2} = rho_w + contribution(salt) + contribution(CO2) */ /* */ /***********************************************************************/ - static Scalar liquidDensity_(Scalar T, - Scalar pl, - Scalar xlH2O, - Scalar xlCO2) + template + static LhsEval liquidDensity_(const LhsEval& T, + const LhsEval& pl, + const LhsEval& xlH2O, + const LhsEval& xlCO2) { Valgrind::CheckDefined(T); Valgrind::CheckDefined(pl); @@ -481,37 +502,40 @@ private: if(T < 273.15) { OPM_THROW(NumericalIssue, - "Liquid density for Brine and CO2 is only " - "defined above 273.15K (is " << T << "K)"); + "Liquid density for Brine and CO2 is only " + "defined above 273.15K (is " << T << "K)"); } if(pl >= 2.5e8) { OPM_THROW(NumericalIssue, - "Liquid density for Brine and CO2 is only " - "defined below 250MPa (is " << pl << "Pa)"); + "Liquid density for Brine and CO2 is only " + "defined below 250MPa (is " << pl << "Pa)"); } - Scalar rho_brine = Brine::liquidDensity(T, pl); - Scalar rho_pure = H2O::liquidDensity(T, pl); - Scalar rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlH2O, xlCO2); - Scalar contribCO2 = rho_lCO2 - rho_pure; + const LhsEval& rho_brine = Brine::liquidDensity(T, pl); + const LhsEval& rho_pure = H2O::liquidDensity(T, pl); + const LhsEval& rho_lCO2 = liquidDensityWaterCO2_(T, pl, xlH2O, xlCO2); + const LhsEval& contribCO2 = rho_lCO2 - rho_pure; return rho_brine + contribCO2; } - static Scalar liquidDensityWaterCO2_(Scalar temperature, - Scalar pl, - Scalar xlH2O, - Scalar xlCO2) + template + static LhsEval liquidDensityWaterCO2_(const LhsEval& temperature, + const LhsEval& pl, + const LhsEval& /*xlH2O*/, + const LhsEval& xlCO2) { - const Scalar M_CO2 = CO2::molarMass(); - const Scalar M_H2O = H2O::molarMass(); + Scalar M_CO2 = CO2::molarMass(); + Scalar M_H2O = H2O::molarMass(); - const Scalar tempC = temperature - 273.15; /* tempC : temperature in °C */ - const Scalar rho_pure = H2O::liquidDensity(temperature, pl); - xlH2O = 1.0 - xlCO2; // xlH2O is available, but in case of a pure gas phase - // the value of M_T for the virtual liquid phase can become very large - const Scalar M_T = M_H2O * xlH2O + M_CO2 * xlCO2; - const Scalar V_phi = + const LhsEval& tempC = temperature - 273.15; /* tempC : temperature in °C */ + const LhsEval& rho_pure = H2O::liquidDensity(temperature, pl); + // calculate the mole fraction of CO2 in the liquid. note that xlH2O is available + // as a function parameter, but in the case of a pure gas phase the value of M_T + // for the virtual liquid phase can become very large + const LhsEval xlH2O = 1.0 - xlCO2; + const LhsEval& M_T = M_H2O * xlH2O + M_CO2 * xlCO2; + const LhsEval& V_phi = (37.51 + tempC*(-9.585e-2 + tempC*(8.74e-4 - @@ -519,41 +543,43 @@ private: return 1/ (xlCO2 * V_phi/M_T + M_H2O * xlH2O / (rho_pure * M_T)); } - static Scalar liquidEnthalpyBrineCO2_(Scalar T, - Scalar p, - Scalar S, - Scalar X_CO2_w) + template + static LhsEval liquidEnthalpyBrineCO2_(const LhsEval& T, + const LhsEval& p, + Scalar S, // salinity + const LhsEval& X_CO2_w) { + typedef MathToolbox LhsToolbox; + /* X_CO2_w : mass fraction of CO2 in brine */ /* same function as enthalpy_brine, only extended by CO2 content */ /*Numerical coefficents from PALLISER*/ - static const Scalar f[] = { + static Scalar f[] = { 2.63500E-1, 7.48368E-6, 1.44611E-6, -3.80860E-10 }; /*Numerical coefficents from MICHAELIDES for the enthalpy of brine*/ - static const Scalar a[4][3] = { + static Scalar a[4][3] = { { 9633.6, -4080.0, +286.49 }, { +166.58, +68.577, -4.6856 }, { -0.90963, -0.36524, +0.249667E-1 }, { +0.17965E-2, +0.71924E-3, -0.4900E-4 } }; - Scalar theta, h_NaCl; - Scalar m, h_ls1, d_h; - Scalar S_lSAT, delta_h; - int i, j; - Scalar delta_hCO2, hg, hw; + LhsEval theta, h_NaCl; + LhsEval h_ls1, d_h; + LhsEval delta_h; + LhsEval delta_hCO2, hg, hw; theta = T - 273.15; - S_lSAT = f[0] + f[1]*theta + f[2]*theta*theta + f[3]*theta*theta*theta; - /*Regularization*/ - if (S>S_lSAT) { + // Regularization + Scalar scalarTheta = LhsToolbox::value(theta); + Scalar S_lSAT = f[0] + scalarTheta*(f[1] + scalarTheta*(f[2] + scalarTheta*f[3])); + if (S > S_lSAT) S = S_lSAT; - } hw = H2O::liquidEnthalpy(T, p) /1E3; /* kJ/kg */ @@ -561,14 +587,14 @@ private: /*U=*/h_NaCl = (3.6710E4*T + 0.5*(6.2770E1)*T*T - ((6.6670E-2)/3)*T*T*T +((2.8000E-5)/4)*(T*T*T*T))/(58.44E3)- 2.045698e+02; /* kJ/kg */ - m = (1E3/58.44)*(S/(1-S)); - i = 0; - j = 0; + Scalar m = 1E3/58.44 * S/(1-S); + int i = 0; + int j = 0; d_h = 0; for (i = 0; i<=3; i++) { for (j=0; j<=2; j++) { - d_h = d_h + a[i][j] * pow(theta, i) * pow(m, j); + d_h = d_h + a[i][j] * LhsToolbox::pow(theta, static_cast(i)) * std::pow(m, j); } } /* heat of dissolution for halite according to Michaelides 1971 */ diff --git a/opm/material/fluidsystems/GasPhase.hpp b/opm/material/fluidsystems/GasPhase.hpp index 4bc818f94..a9375494b 100644 --- a/opm/material/fluidsystems/GasPhase.hpp +++ b/opm/material/fluidsystems/GasPhase.hpp @@ -66,25 +66,25 @@ public: * \brief The mass in [kg] of one mole of the component. */ static Scalar molarMass() - { return Component::molarMass(); } + { return Component::molarMass(); } /*! * \brief Returns the critical temperature of the component */ static Scalar criticalTemperature() - { return Component::criticalTemperature(); } + { return Component::criticalTemperature(); } /*! * \brief Returns the critical pressure of the component */ static Scalar criticalPressure() - { return Component::criticalPressure(); } + { return Component::criticalPressure(); } /*! * \brief Returns the temperature at the component's triple point. */ static Scalar tripleTemperature() - { return Component::tripleTemperature(); } + { return Component::tripleTemperature(); } /*! * \brief Returns the pressure at the component's triple point. @@ -98,7 +98,8 @@ public: * * \copydetails Doxygen::temperatureParam */ - static Scalar vaporPressure(Scalar temperature) + template + static Evaluation vaporPressure(const Evaluation& temperature) { return Component::vaporPressure(temperature); } /*! @@ -106,8 +107,9 @@ public: * * \copydetails Doxygen::TpParams */ - static Scalar density(Scalar temperature, Scalar pressure) - { return Component::gasDensity(temperature, pressure); } + template + static Evaluation density(const Evaluation& temperature, const Evaluation& pressure) + { return Component::gasDensity(temperature, pressure); } /*! * \brief The pressure [Pa] of the component at a given density and temperature. @@ -115,23 +117,26 @@ public: * \param temperature The temperature of interest [K] * \param density The density of interest [kg / m^3] */ - static Scalar pressure(Scalar temperature, Scalar density) - { return Component::gasPressure(temperature, density); } + template + static Evaluation pressure(const Evaluation& temperature, const Evaluation& density) + { return Component::gasPressure(temperature, density); } /*! * \brief Specific enthalpy [J/kg] the pure component as a gas. * * \copydetails Doxygen::TpParams */ - static const Scalar enthalpy(Scalar temperature, Scalar pressure) - { return Component::gasEnthalpy(temperature, pressure); } + template + static Evaluation enthalpy(const Evaluation& temperature, const Evaluation& pressure) + { return Component::gasEnthalpy(temperature, pressure); } /*! * \brief Specific internal energy [J/kg] the pure component as a gas. * * \copydetails Doxygen::TpParams */ - static const Scalar internalEnergy(Scalar temperature, Scalar pressure) + template + static Evaluation internalEnergy(const Evaluation& temperature, const Evaluation& pressure) { return Component::gasInternalEnergy(temperature, pressure); } /*! @@ -139,15 +144,17 @@ public: * * \copydetails Doxygen::TpParams */ - static Scalar viscosity(Scalar temperature, Scalar pressure) - { return Component::gasViscosity(temperature, pressure); } + template + static Evaluation viscosity(const Evaluation& temperature, const Evaluation& pressure) + { return Component::gasViscosity(temperature, pressure); } /*! * \brief Thermal conductivity of the fluid [W/(m^2 K/m)]. * * \copydetails Doxygen::TpParams */ - static Scalar thermalConductivity(Scalar temperature, Scalar pressure) + template + static Evaluation thermalConductivity(const Evaluation& temperature, const Evaluation& pressure) { return Component::gasThermalConductivity(temperature, pressure); } /*! @@ -155,7 +162,8 @@ public: * * \copydetails Doxygen::TpParams */ - static Scalar heatCapacity(Scalar temperature, Scalar pressure) + template + static Evaluation heatCapacity(const Evaluation& temperature, const Evaluation& pressure) { return Component::gasHeatCapacity(temperature, pressure); } }; } // namespace Opm diff --git a/opm/material/fluidsystems/H2OAirFluidSystem.hpp b/opm/material/fluidsystems/H2OAirFluidSystem.hpp index 578ca12c4..5781d3e7f 100644 --- a/opm/material/fluidsystems/H2OAirFluidSystem.hpp +++ b/opm/material/fluidsystems/H2OAirFluidSystem.hpp @@ -249,17 +249,20 @@ public: } //! \copydoc BaseFluidSystem::density - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + typedef Opm::MathToolbox LhsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p; + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + LhsEval p; if (isCompressible(phaseIdx)) - p = fluidState.pressure(phaseIdx); + p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); else { // random value which will hopefully cause things to blow // up if it is used in a calculation! @@ -268,9 +271,9 @@ public: } - Scalar sumMoleFrac = 0; + LhsEval sumMoleFrac = 0; for (int compIdx = 0; compIdx < numComponents; ++compIdx) - sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx); + sumMoleFrac += FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, compIdx)); if (phaseIdx == liquidPhaseIdx) { @@ -280,15 +283,13 @@ public: else { // See: Ochs 2008 (2.6) - Scalar rholH2O = H2O::liquidDensity(T, p); - Scalar clH2O = rholH2O/H2O::molarMass(); + const LhsEval& rholH2O = H2O::liquidDensity(T, p); + const LhsEval& clH2O = rholH2O/H2O::molarMass(); - return - clH2O - * (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx) - + - Air::molarMass()*fluidState.moleFraction(liquidPhaseIdx, AirIdx)) - / sumMoleFrac; + const auto& xlH2O = FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, H2OIdx)); + const auto& xlAir = FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, AirIdx)); + + return clH2O*(H2O::molarMass()*xlH2O + Air::molarMass()*xlAir)/sumMoleFrac; } } else if (phaseIdx == gasPhaseIdx) @@ -297,34 +298,35 @@ public: // for the gas phase assume an ideal gas return IdealGas::molarDensity(T, p) - * fluidState.averageMolarMass(gasPhaseIdx) - / std::max(1e-5, sumMoleFrac); + * FsToolbox::template toLhs(fluidState.averageMolarMass(gasPhaseIdx)) + / LhsToolbox::max(1e-5, sumMoleFrac); - Scalar partialPressureH2O = - fluidState.moleFraction(gasPhaseIdx, H2OIdx) * - fluidState.pressure(gasPhaseIdx); + LhsEval partialPressureH2O = + FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx)) + *FsToolbox::template toLhs(fluidState.pressure(gasPhaseIdx)); - Scalar partialPressureAir = - fluidState.moleFraction(gasPhaseIdx, AirIdx) * - fluidState.pressure(gasPhaseIdx); + LhsEval partialPressureAir = + FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, AirIdx)) + *FsToolbox::template toLhs(fluidState.pressure(gasPhaseIdx)); - return - H2O::gasDensity(T, partialPressureH2O) + - Air::gasDensity(T, partialPressureAir); + return H2O::gasDensity(T, partialPressureH2O) + Air::gasDensity(T, partialPressureAir); } OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, + template + static LhsEval viscosity(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { + typedef Opm::MathToolbox LhsToolbox; + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == liquidPhaseIdx) { @@ -348,10 +350,9 @@ public: * */ - Scalar muResult = 0; - const Scalar mu[numComponents] = { - H2O::gasViscosity(T, - H2O::vaporPressure(T)), + LhsEval muResult = 0; + const LhsEval mu[numComponents] = { + H2O::gasViscosity(T, H2O::vaporPressure(T)), Air::gasViscosity(T, p) }; @@ -363,17 +364,20 @@ public: for (int i = 0; i < numComponents; ++i) { - Scalar divisor = 0; + LhsEval divisor = 0; for (int j = 0; j < numComponents; ++j) { - Scalar phiIJ = 1 + std::sqrt(mu[i]/mu[j]) * // 1 + (mu[i]/mu[j]^1/2 + LhsEval phiIJ = + 1 + + LhsToolbox::sqrt(mu[i]/mu[j]) * // 1 + (mu[i]/mu[j]^1/2 std::pow(M[j]/M[i], 1./4.0); // (M[i]/M[j])^1/4 phiIJ *= phiIJ; phiIJ /= std::sqrt(8*(1 + M[i]/M[j])); - divisor += fluidState.moleFraction(phaseIdx, j)*phiIJ; + divisor += FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, j))*phiIJ; } - muResult += fluidState.moleFraction(phaseIdx, i)*mu[i] / divisor; + const auto& xAlphaI = FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, i)); + muResult += xAlphaI*mu[i]/divisor; } return muResult; } @@ -382,17 +386,19 @@ public: } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx, int compIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == liquidPhaseIdx) { if (compIdx == H2OIdx) @@ -406,14 +412,16 @@ public: } //! \copydoc BaseFluidSystem::diffusionCoefficient - template - static Scalar binaryDiffusionCoefficient(const FluidState &fluidState, + template + static LhsEval binaryDiffusionCoefficient(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx, int compIdx) { - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == liquidPhaseIdx) return BinaryCoeff::H2O_Air::liquidDiffCoeff(T, p); @@ -423,13 +431,15 @@ public: } //! \copydoc BaseFluidSystem::enthalpy - template - static Scalar enthalpy(const FluidState &fluidState, + template + static LhsEval enthalpy(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); Valgrind::CheckDefined(T); Valgrind::CheckDefined(p); @@ -441,48 +451,52 @@ public: else if (phaseIdx == gasPhaseIdx) { - Scalar result = 0.0; + LhsEval result = 0.0; result += H2O::gasEnthalpy(T, p) * - fluidState.massFraction(gasPhaseIdx, H2OIdx); + FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, H2OIdx)); result += Air::gasEnthalpy(T, p) * - fluidState.massFraction(gasPhaseIdx, AirIdx); + FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, AirIdx)); return result; } OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx); } //! \copydoc BaseFluidSystem::thermalConductivity - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx) ; - Scalar pressure = fluidState.pressure(phaseIdx); + const LhsEval& temperature = + FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& pressure = + FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); - if (phaseIdx == liquidPhaseIdx){// liquid phase + if (phaseIdx == liquidPhaseIdx) return H2O::liquidThermalConductivity(temperature, pressure); - } - else{// gas phase - Scalar lambdaDryAir = Air::gasThermalConductivity(temperature, pressure); + else { // gas phase + const LhsEval& lambdaDryAir = Air::gasThermalConductivity(temperature, pressure); if (useComplexRelations){ - Scalar xAir = fluidState.moleFraction(phaseIdx, AirIdx); - Scalar xH2O = fluidState.moleFraction(phaseIdx, H2OIdx); - Scalar lambdaAir = xAir * lambdaDryAir; + const LhsEval& xAir = + FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, AirIdx)); + const LhsEval& xH2O = + FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, H2OIdx)); + LhsEval lambdaAir = xAir*lambdaDryAir; // Assuming Raoult's, Daltons law and ideal gas // in order to obtain the partial density of water in the air phase - Scalar partialPressure = pressure * xH2O; + LhsEval partialPressure = pressure*xH2O; - Scalar lambdaH2O = - xH2O - * H2O::gasThermalConductivity(temperature, partialPressure); + LhsEval lambdaH2O = + xH2O*H2O::gasThermalConductivity(temperature, partialPressure); return lambdaAir + lambdaH2O; } else diff --git a/opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp b/opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp index 2128308ec..d4dbff112 100644 --- a/opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp +++ b/opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp @@ -194,44 +194,46 @@ public: } //! \copydoc BaseFluidSystem::density - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { - Scalar T = fluidState.temperature(phaseIdx) ; + typedef Opm::MathToolbox FsToolbox; + + const LhsEval& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); if (phaseIdx == waterPhaseIdx) { // See: Ochs 2008 - Scalar p = H2O::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100; - Scalar rholH2O = H2O::liquidDensity(T, p); - Scalar clH2O = rholH2O/H2O::molarMass(); + const LhsEval& p = + H2O::liquidIsCompressible() + ? FsToolbox::template toLhs(fluidState.pressure(phaseIdx)) + : 1e100; + + const LhsEval& rholH2O = H2O::liquidDensity(T, p); + const LhsEval& clH2O = rholH2O/H2O::molarMass(); // this assumes each dissolved molecule displaces exactly one // water molecule in the liquid return - clH2O*(H2O::molarMass()*fluidState.moleFraction(waterPhaseIdx, H2OIdx) - + - Air::molarMass()*fluidState.moleFraction(waterPhaseIdx, airIdx) - + - NAPL::molarMass()*fluidState.moleFraction(waterPhaseIdx, NAPLIdx)); + clH2O*(H2O::molarMass()*FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, H2OIdx)) + + Air::molarMass()*FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, airIdx)) + + NAPL::molarMass()*FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, NAPLIdx))); } else if (phaseIdx == naplPhaseIdx) { // assume pure NAPL for the NAPL phase - Scalar p = NAPL::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100; + const LhsEval& p = + NAPL::liquidIsCompressible() + ? FsToolbox::template toLhs(fluidState.pressure(phaseIdx)) + : 1e100; return NAPL::liquidDensity(T, p); } assert (phaseIdx == gasPhaseIdx); - Scalar pH2O = - fluidState.moleFraction(gasPhaseIdx, H2OIdx) * - fluidState.pressure(gasPhaseIdx); - Scalar pAir = - fluidState.moleFraction(gasPhaseIdx, airIdx) * - fluidState.pressure(gasPhaseIdx); - Scalar pNAPL = - fluidState.moleFraction(gasPhaseIdx, NAPLIdx) * - fluidState.pressure(gasPhaseIdx); + const LhsEval& pg = FsToolbox::template toLhs(fluidState.pressure(gasPhaseIdx)); + const LhsEval& pH2O = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*pg; + const LhsEval& pAir = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, airIdx))*pg; + const LhsEval& pNAPL = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*pg; return H2O::gasDensity(T, pH2O) + Air::gasDensity(T, pAir) + @@ -239,13 +241,15 @@ public: } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const LhsEval& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == waterPhaseIdx) { // assume pure water viscosity @@ -271,8 +275,7 @@ public: * divisions * -- compare e.g. with Promo Class p. 32/33 */ - Scalar muResult; - const Scalar mu[numComponents] = { + const LhsEval mu[numComponents] = { H2O::gasViscosity(T, H2O::vaporPressure(T)), Air::gasViscosity(T, p), NAPL::gasViscosity(T, NAPL::vaporPressure(T)) @@ -284,66 +287,61 @@ public: NAPL::molarMass() }; - Scalar muAW = mu[airIdx]*fluidState.moleFraction(gasPhaseIdx, airIdx) - + mu[H2OIdx]*fluidState.moleFraction(gasPhaseIdx, H2OIdx) - / (fluidState.moleFraction(gasPhaseIdx, airIdx) - + fluidState.moleFraction(gasPhaseIdx, H2OIdx)); - Scalar xAW = fluidState.moleFraction(gasPhaseIdx, airIdx) - + fluidState.moleFraction(gasPhaseIdx, H2OIdx); - - Scalar MAW = (fluidState.moleFraction(gasPhaseIdx, airIdx)*Air::molarMass() - + fluidState.moleFraction(gasPhaseIdx, H2OIdx)*H2O::molarMass()) - / xAW; + const LhsEval& xgAir = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, airIdx)); + const LhsEval& xgH2O = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx)); + const LhsEval& xgNapl = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, NAPLIdx)); + const LhsEval& xgAW = xgAir + xgH2O; + const LhsEval& muAW = (mu[airIdx]*xgAir + mu[H2OIdx]*xgH2O)/xgAW; + const LhsEval& MAW = (xgAir*Air::molarMass() + xgH2O*H2O::molarMass())/xgAW; Scalar phiCAW = 0.3; // simplification for this particular system /* actually like this * Scalar phiCAW = std::pow(1.+std::sqrt(mu[NAPLIdx]/muAW)*std::pow(MAW/M[NAPLIdx],0.25),2) * / std::sqrt(8.*(1.+M[NAPLIdx]/MAW)); */ - Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW); + const LhsEval& phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW); - muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gasPhaseIdx, NAPLIdx)*phiAWC) - + (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) * mu[NAPLIdx]) - / (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) + xAW*phiCAW); - return muResult; + return (xgAW*muAW)/(xgAW + xgNapl*phiAWC) + (xgNapl*mu[NAPLIdx])/(xgNapl + xgAW*phiCAW); } //! \copydoc BaseFluidSystem::diffusionCoefficient - template - static Scalar diffusionCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval diffusionCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { return 0; #if 0 - Scalar T = fluidState.temperature(phaseIdx) ; - Scalar p = fluidState.pressure(phaseIdx); - Scalar diffCont; + typedef Opm::MathToolbox FsToolbox; + + const LhsEval& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + LhsEval diffCont; if (phaseIdx==gasPhaseIdx) { - Scalar diffAC = Opm::BinaryCoeff::Air_Mesitylene::gasDiffCoeff(T, p); - Scalar diffWC = Opm::BinaryCoeff::H2O_Mesitylene::gasDiffCoeff(T, p); - Scalar diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p); + const LhsEval& diffAC = Opm::BinaryCoeff::Air_Mesitylene::gasDiffCoeff(T, p); + const LhsEval& diffWC = Opm::BinaryCoeff::H2O_Mesitylene::gasDiffCoeff(T, p); + const LhsEval& diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p); - const Scalar xga = fluidState.moleFraction(gasPhaseIdx, airIdx); - const Scalar xgw = fluidState.moleFraction(gasPhaseIdx, H2OIdx); - const Scalar xgc = fluidState.moleFraction(gasPhaseIdx, NAPLIdx); + const LhsEval& xga = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, airIdx)); + const LhsEval& xgw = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx)); + const LhsEval& xgc = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, NAPLIdx)); if (compIdx==NAPLIdx) return (1 - xgw)/(xga/diffAW + xgc/diffWC); else if (compIdx==H2OIdx) return (1 - xgc)/(xgw/diffWC + xga/diffAC); else if (compIdx==airIdx) OPM_THROW(std::logic_error, - "Diffusivity of air in the gas phase " - "is constraint by sum of diffusive fluxes = 0 !\n"); + "Diffusivity of air in the gas phase " + "is constraint by sum of diffusive fluxes = 0 !\n"); } else if (phaseIdx==waterPhaseIdx){ - Scalar diffACl = 1.e-9; // BinaryCoeff::Air_Mesitylene::liquidDiffCoeff(temperature, pressure); - Scalar diffWCl = 1.e-9; // BinaryCoeff::H2O_Mesitylene::liquidDiffCoeff(temperature, pressure); - Scalar diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure); + const LhsEval& diffACl = 1.e-9; // BinaryCoeff::Air_Mesitylene::liquidDiffCoeff(temperature, pressure); + const LhsEval& diffWCl = 1.e-9; // BinaryCoeff::H2O_Mesitylene::liquidDiffCoeff(temperature, pressure); + const LhsEval& diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure); - Scalar xwa = fluidState.moleFraction(waterPhaseIdx, airIdx); - Scalar xww = fluidState.moleFraction(waterPhaseIdx, H2OIdx); - Scalar xwc = fluidState.moleFraction(waterPhaseIdx, NAPLIdx); + const LhsEval& xwa = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, airIdx)); + const LhsEval& xww = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, H2OIdx)); + const LhsEval& xwc = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, NAPLIdx)); switch (compIdx) { case NAPLIdx: @@ -354,31 +352,33 @@ public: return diffCont; case H2OIdx: OPM_THROW(std::logic_error, - "Diffusivity of water in the water phase " - "is constraint by sum of diffusive fluxes = 0 !\n"); + "Diffusivity of water in the water phase " + "is constraint by sum of diffusive fluxes = 0 !\n"); }; } else if (phaseIdx==naplPhaseIdx) { OPM_THROW(std::logic_error, - "Diffusion coefficients of " - "substances in liquid phase are undefined!\n"); + "Diffusion coefficients of " + "substances in liquid phase are undefined!\n"); } return 0; #endif } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const LhsEval& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); Valgrind::CheckDefined(T); Valgrind::CheckDefined(p); @@ -397,7 +397,7 @@ public: // other components, i.e. the fugacity cofficient is much // smaller. else if (phaseIdx == naplPhaseIdx) { - Scalar phiNapl = NAPL::vaporPressure(T)/p; + const LhsEval& phiNapl = NAPL::vaporPressure(T)/p; if (compIdx == NAPLIdx) return phiNapl; else if (compIdx == airIdx) @@ -415,13 +415,15 @@ public: //! \copydoc BaseFluidSystem::enthalpy - template - static Scalar enthalpy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { - Scalar T = fluidState.temperature(phaseIdx) ; - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const LhsEval& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == waterPhaseIdx) { return H2O::liquidEnthalpy(T, p); @@ -431,10 +433,10 @@ public: } else if (phaseIdx == gasPhaseIdx) { // gas phase enthalpy depends strongly on composition - Scalar result = 0; - result += H2O::gasEnthalpy(T, p) * fluidState.massFraction(gasPhaseIdx, H2OIdx); - result += NAPL::gasEnthalpy(T, p) * fluidState.massFraction(gasPhaseIdx, airIdx); - result += Air::gasEnthalpy(T, p) * fluidState.massFraction(gasPhaseIdx, NAPLIdx); + LhsEval result = 0; + result += H2O::gasEnthalpy(T, p) * FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, H2OIdx)); + result += NAPL::gasEnthalpy(T, p) * FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, airIdx)); + result += Air::gasEnthalpy(T, p) * FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, NAPLIdx)); return result; } @@ -442,22 +444,26 @@ public: } //! \copydoc BaseFluidSystem::thermalConductivity - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar T = fluidState.temperature(phaseIdx) ; - Scalar p = fluidState.pressure(phaseIdx); - if (phaseIdx == waterPhaseIdx){ // water phase + const LhsEval& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + return H2O::liquidThermalConductivity(T, p); } else if (phaseIdx == gasPhaseIdx) { // gas phase - Scalar lambdaDryAir = Air::gasThermalConductivity(T, p); - return lambdaDryAir; + const LhsEval& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const LhsEval& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + + return Air::gasThermalConductivity(T, p); } assert(phaseIdx == naplPhaseIdx); diff --git a/opm/material/fluidsystems/H2OAirXyleneFluidSystem.hpp b/opm/material/fluidsystems/H2OAirXyleneFluidSystem.hpp index 2b4923329..67234b98b 100644 --- a/opm/material/fluidsystems/H2OAirXyleneFluidSystem.hpp +++ b/opm/material/fluidsystems/H2OAirXyleneFluidSystem.hpp @@ -161,62 +161,66 @@ public: } //! \copydoc BaseFluidSystem::density - template - static Scalar density(const FluidState &fluidState, + template + static LhsEval density(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + if (phaseIdx == waterPhaseIdx) { + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + // See: Ochs 2008 // \todo: proper citation - Scalar rholH2O = H2O::liquidDensity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); - Scalar clH2O = rholH2O/H2O::molarMass(); + const LhsEval& rholH2O = H2O::liquidDensity(T, p); + const LhsEval& clH2O = rholH2O/H2O::molarMass(); + const auto& xwH2O = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, H2OIdx)); + const auto& xwAir = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, airIdx)); + const auto& xwNapl = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, NAPLIdx)); // this assumes each dissolved molecule displaces exactly one // water molecule in the liquid - return - clH2O*(H2O::molarMass()*fluidState.moleFraction(waterPhaseIdx, H2OIdx) - + - Air::molarMass()*fluidState.moleFraction(waterPhaseIdx, airIdx) - + - NAPL::molarMass()*fluidState.moleFraction(waterPhaseIdx, NAPLIdx)); + return clH2O*(H2O::molarMass()*xwH2O + Air::molarMass()*xwAir + NAPL::molarMass()*xwNapl); } else if (phaseIdx == naplPhaseIdx) { // assume pure NAPL for the NAPL phase - Scalar pressure = NAPL::liquidIsCompressible()?fluidState.pressure(phaseIdx):1e100; - return NAPL::liquidDensity(fluidState.temperature(phaseIdx), pressure); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + return NAPL::liquidDensity(T, LhsEval(1e100)); } assert (phaseIdx == gasPhaseIdx); - Scalar pH2O = - fluidState.moleFraction(gasPhaseIdx, H2OIdx) * - fluidState.pressure(gasPhaseIdx); - Scalar pAir = - fluidState.moleFraction(gasPhaseIdx, airIdx) * - fluidState.pressure(gasPhaseIdx); - Scalar pNAPL = - fluidState.moleFraction(gasPhaseIdx, NAPLIdx) * - fluidState.pressure(gasPhaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + + const LhsEval& pH2O = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*p; + const LhsEval& pAir = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, airIdx))*p; + const LhsEval& pNAPL = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*p; return - H2O::gasDensity(fluidState.temperature(phaseIdx), pH2O) + - Air::gasDensity(fluidState.temperature(phaseIdx), pAir) + - NAPL::gasDensity(fluidState.temperature(phaseIdx), pNAPL); + H2O::gasDensity(T, pH2O) + + Air::gasDensity(T, pAir) + + NAPL::gasDensity(T, pNAPL); } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, + template + static LhsEval viscosity(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + if (phaseIdx == waterPhaseIdx) { // assume pure water viscosity - return H2O::liquidViscosity(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx)); + return H2O::liquidViscosity(T, p); } else if (phaseIdx == naplPhaseIdx) { // assume pure NAPL viscosity - return NAPL::liquidViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + return NAPL::liquidViscosity(T, p); } assert (phaseIdx == gasPhaseIdx); @@ -232,11 +236,10 @@ public: * divisions * -- compare e.g. with Promo Class p. 32/33 */ - Scalar muResult; - const Scalar mu[numComponents] = { - H2O::gasViscosity(fluidState.temperature(phaseIdx), H2O::vaporPressure(fluidState.temperature(phaseIdx))), - Air::simpleGasViscosity(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)), - NAPL::gasViscosity(fluidState.temperature(phaseIdx), NAPL::vaporPressure(fluidState.temperature(phaseIdx))) + const LhsEval mu[numComponents] = { + H2O::gasViscosity(T, H2O::vaporPressure(T)), + Air::simpleGasViscosity(T, p), + NAPL::gasViscosity(T, NAPL::vaporPressure(T)) }; // molar masses const Scalar M[numComponents] = { @@ -245,50 +248,45 @@ public: NAPL::molarMass() }; - Scalar muAW = mu[airIdx]*fluidState.moleFraction(gasPhaseIdx, airIdx) - + mu[H2OIdx]*fluidState.moleFraction(gasPhaseIdx, H2OIdx) - / (fluidState.moleFraction(gasPhaseIdx, airIdx) - + fluidState.moleFraction(gasPhaseIdx, H2OIdx)); - Scalar xAW = fluidState.moleFraction(gasPhaseIdx, airIdx) - + fluidState.moleFraction(gasPhaseIdx, H2OIdx); + const auto& xgAir = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, airIdx)); + const auto& xgH2O = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx)); + const auto& xgNapl = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, NAPLIdx)); - Scalar MAW = (fluidState.moleFraction(gasPhaseIdx, airIdx)*Air::molarMass() - + fluidState.moleFraction(gasPhaseIdx, H2OIdx)*H2O::molarMass()) - / xAW; - - /* TODO, please check phiCAW for the Xylene case here */ + const LhsEval& xgAW = xgAir + xgH2O; + const LhsEval& muAW = (mu[airIdx]*xgAir + mu[H2OIdx]*xgH2O)/ xgAW; + const LhsEval& MAW = (xgAir*Air::molarMass() + xgH2O*H2O::molarMass())/xgAW; Scalar phiCAW = 0.3; // simplification for this particular system /* actually like this * Scalar phiCAW = std::pow(1.+std::sqrt(mu[NAPLIdx]/muAW)*std::pow(MAW/M[NAPLIdx],0.25),2) * / std::sqrt(8.*(1.+M[NAPLIdx]/MAW)); */ - Scalar phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW); + const LhsEval& phiAWC = phiCAW * muAW*M[NAPLIdx]/(mu[NAPLIdx]*MAW); - muResult = (xAW*muAW)/(xAW+fluidState.moleFraction(gasPhaseIdx, NAPLIdx)*phiAWC) - + (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) * mu[NAPLIdx]) - / (fluidState.moleFraction(gasPhaseIdx, NAPLIdx) + xAW*phiCAW); - return muResult; + return (xgAW*muAW)/(xgAW+xgNapl*phiAWC) + (xgNapl*mu[NAPLIdx])/(xgNapl + xgAW*phiCAW); } //! \copydoc BaseFluidSystem::diffusionCoefficient - template - static Scalar diffusionCoefficient(const FluidState &fluidState, + template + static LhsEval diffusionCoefficient(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx, int compIdx) { - Scalar diffCont; + typedef Opm::MathToolbox FsToolbox; if (phaseIdx==gasPhaseIdx) { - Scalar diffAC = Opm::BinaryCoeff::Air_Xylene::gasDiffCoeff(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); - Scalar diffWC = Opm::BinaryCoeff::H2O_Xylene::gasDiffCoeff(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); - Scalar diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); - const Scalar xga = fluidState.moleFraction(gasPhaseIdx, airIdx); - const Scalar xgw = fluidState.moleFraction(gasPhaseIdx, H2OIdx); - const Scalar xgc = fluidState.moleFraction(gasPhaseIdx, NAPLIdx); + const LhsEval& diffAC = Opm::BinaryCoeff::Air_Xylene::gasDiffCoeff(T, p); + const LhsEval& diffWC = Opm::BinaryCoeff::H2O_Xylene::gasDiffCoeff(T, p); + const LhsEval& diffAW = Opm::BinaryCoeff::H2O_Air::gasDiffCoeff(T, p); + + const LhsEval& xga = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, airIdx)); + const LhsEval& xgw = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx)); + const LhsEval& xgc = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, NAPLIdx)); if (compIdx==NAPLIdx) return (1.- xgw)/(xga/diffAW + xgc/diffWC); else if (compIdx==H2OIdx) return (1.- xgc)/(xgw/diffWC + xga/diffAC); @@ -300,17 +298,15 @@ public: Scalar diffWCl = 1.e-9; // BinaryCoeff::H2O_Xylene::liquidDiffCoeff(temperature, pressure); Scalar diffAWl = 1.e-9; // BinaryCoeff::H2O_Air::liquidDiffCoeff(temperature, pressure); - Scalar xwa = fluidState.moleFraction(waterPhaseIdx, airIdx); - Scalar xww = fluidState.moleFraction(waterPhaseIdx, H2OIdx); - Scalar xwc = fluidState.moleFraction(waterPhaseIdx, NAPLIdx); + const LhsEval& xwa = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, airIdx)); + const LhsEval& xww = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, H2OIdx)); + const LhsEval& xwc = FsToolbox::template toLhs(fluidState.moleFraction(waterPhaseIdx, NAPLIdx)); switch (compIdx) { case NAPLIdx: - diffCont = (1.- xww)/(xwa/diffAWl + xwc/diffWCl); - return diffCont; + return (1.- xww)/(xwa/diffAWl + xwc/diffWCl); case airIdx: - diffCont = (1.- xwc)/(xww/diffWCl + xwa/diffACl); - return diffCont; + return (1.- xwc)/(xww/diffWCl + xwa/diffACl); case H2OIdx: OPM_THROW(std::logic_error, "Diffusivity of water in the water phase " @@ -326,17 +322,19 @@ public: } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx, int compIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == waterPhaseIdx) { if (compIdx == H2OIdx) @@ -353,7 +351,7 @@ public: // other components, i.e. the fugacity cofficient is much // smaller. if (phaseIdx == naplPhaseIdx) { - Scalar phiNapl = NAPL::vaporPressure(T)/p; + const LhsEval& phiNapl = NAPL::vaporPressure(T)/p; if (compIdx == NAPLIdx) return phiNapl; else if (compIdx == airIdx) @@ -369,28 +367,31 @@ public: } //! \copydoc BaseFluidSystem::enthalpy - template - static Scalar enthalpy(const FluidState &fluidState, + template + static LhsEval enthalpy(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + if (phaseIdx == waterPhaseIdx) { - return H2O::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + return H2O::liquidEnthalpy(T, p); } else if (phaseIdx == naplPhaseIdx) { - return NAPL::liquidEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); + return NAPL::liquidEnthalpy(T, p); } else if (phaseIdx == gasPhaseIdx) { // gas phase enthalpy depends strongly on composition - Scalar hgc = NAPL::gasEnthalpy(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx)); - Scalar hgw = H2O::gasEnthalpy(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx)); - Scalar hga = Air::gasEnthalpy(fluidState.temperature(phaseIdx), fluidState.pressure(phaseIdx)); // pressure is only a dummy here (not dependent on pressure, just temperature) + const LhsEval& hgc = NAPL::gasEnthalpy(T, p); + const LhsEval& hgw = H2O::gasEnthalpy(T, p); + const LhsEval& hga = Air::gasEnthalpy(T, p); - Scalar result = 0; - result += hgw * fluidState.massFraction(gasPhaseIdx, H2OIdx); - result += hga * fluidState.massFraction(gasPhaseIdx, airIdx); - result += hgc * fluidState.massFraction(gasPhaseIdx, NAPLIdx); + LhsEval result = 0; + result += hgw * FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, H2OIdx)); + result += hga * FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, airIdx)); + result += hgc * FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, NAPLIdx)); return result; } @@ -398,28 +399,32 @@ public: } private: - static Scalar waterPhaseDensity_(Scalar T, Scalar pw, Scalar xww, Scalar xwa, Scalar xwc) + template + static LhsEval waterPhaseDensity_(const LhsEval& T, + const LhsEval& pw, + const LhsEval& xww, + const LhsEval& xwa, + const LhsEval& xwc) { - Scalar rholH2O = H2O::liquidDensity(T, pw); - Scalar clH2O = rholH2O/H2O::molarMass(); + const LhsEval& rholH2O = H2O::liquidDensity(T, pw); + const LhsEval& clH2O = rholH2O/H2O::molarMass(); // this assumes each dissolved molecule displaces exactly one // water molecule in the liquid - return - clH2O*(xww*H2O::molarMass() + xwa*Air::molarMass() + xwc*NAPL::molarMass()); + return clH2O*(xww*H2O::molarMass() + xwa*Air::molarMass() + xwc*NAPL::molarMass()); } - static Scalar gasPhaseDensity_(Scalar T, Scalar pg, Scalar xgw, Scalar xga, Scalar xgc) - { - return H2O::gasDensity(T, pg*xgw) + Air::gasDensity(T, pg*xga) + NAPL::gasDensity(T, pg*xgc); - } - - static Scalar NAPLPhaseDensity_(Scalar T, Scalar pn) - { - return - NAPL::liquidDensity(T, pn); - } + template + static LhsEval gasPhaseDensity_(const LhsEval& T, + const LhsEval& pg, + const LhsEval& xgw, + const LhsEval& xga, + const LhsEval& xgc) + { return H2O::gasDensity(T, pg*xgw) + Air::gasDensity(T, pg*xga) + NAPL::gasDensity(T, pg*xgc); } + template + static LhsEval NAPLPhaseDensity_(const LhsEval& T, const LhsEval& pn) + { return NAPL::liquidDensity(T, pn); } }; } // namespace FluidSystems } // namespace Opm diff --git a/opm/material/fluidsystems/H2ON2FluidSystem.hpp b/opm/material/fluidsystems/H2ON2FluidSystem.hpp index d7d69c1a9..d48ccdca9 100644 --- a/opm/material/fluidsystems/H2ON2FluidSystem.hpp +++ b/opm/material/fluidsystems/H2ON2FluidSystem.hpp @@ -256,26 +256,26 @@ public: /*! * \copydoc BaseFluidSystem::density * - * If useComplexRelations == true, we apply - * Formula (2.6) from S.O.Ochs: - * "Development of a multiphase multicomponent - * model for PEMFC - Technical report: IRTG-NUPUS", + * If useComplexRelations == true, we apply Formula (2.6) from S.O.Ochs: "Development + * of a multiphase multicomponent model for PEMFC - Technical report: IRTG-NUPUS", * University of Stuttgart, 2008 */ - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox LhsToolbox; + typedef Opm::MathToolbox FsToolbox; - Scalar sumMoleFrac = 0; + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + + LhsEval sumMoleFrac = 0; for (int compIdx = 0; compIdx < numComponents; ++compIdx) - sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx); - Valgrind::CheckDefined(sumMoleFrac); + sumMoleFrac += FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, compIdx)); // liquid phase if (phaseIdx == liquidPhaseIdx) { @@ -285,103 +285,112 @@ public: else { // See: Ochs 2008 - Scalar rholH2O = H2O::liquidDensity(T, p); - Scalar clH2O = rholH2O/H2O::molarMass(); + const auto& rholH2O = H2O::liquidDensity(T, p); + const auto& clH2O = rholH2O/H2O::molarMass(); + + const auto& xlH2O = FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, H2OIdx)); + const auto& xlN2 = FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, N2Idx)); // this assumes each nitrogen molecule displaces exactly one // water molecule in the liquid - return - clH2O - * (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx) - + - N2::molarMass()*fluidState.moleFraction(liquidPhaseIdx, N2Idx)) - / sumMoleFrac; + return clH2O*(H2O::molarMass()*xlH2O + N2::molarMass()*xlN2)/sumMoleFrac; } } // gas phase + assert(phaseIdx == gasPhaseIdx); + if (!useComplexRelations) // for the gas phase assume an ideal gas return IdealGas::molarDensity(T, p) - * fluidState.averageMolarMass(gasPhaseIdx) - / std::max(1e-5, sumMoleFrac); + * FsToolbox::template toLhs(fluidState.averageMolarMass(gasPhaseIdx)) + / LhsToolbox::max(1e-5, sumMoleFrac); // assume ideal mixture: steam and nitrogen don't "see" each // other - Scalar rho_gH2O = H2O::gasDensity(T, p*fluidState.moleFraction(gasPhaseIdx, H2OIdx)); - Scalar rho_gN2 = N2::gasDensity(T, p*fluidState.moleFraction(gasPhaseIdx, N2Idx)); - return (rho_gH2O + rho_gN2) / std::max(1e-5, sumMoleFrac); + const auto& xgH2O = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, H2OIdx)); + const auto& xgN2 = FsToolbox::template toLhs(fluidState.moleFraction(gasPhaseIdx, N2Idx)); + const auto& rho_gH2O = H2O::gasDensity(T, p*xgH2O); + const auto& rho_gN2 = N2::gasDensity(T, p*xgN2); + return (rho_gH2O + rho_gN2)/LhsToolbox::max(1e-5, sumMoleFrac); } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox LhsToolbox; + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); // liquid phase - if (phaseIdx == liquidPhaseIdx) { + if (phaseIdx == liquidPhaseIdx) // assume pure water for the liquid phase return H2O::liquidViscosity(T, p); - } // gas phase + assert(phaseIdx == gasPhaseIdx); + if (!useComplexRelations) - { // assume pure nitrogen for the gas phase return N2::gasViscosity(T, p); - } - else - { + else { /* Wilke method. See: * * See: R. Reid, et al.: The Properties of Gases and Liquids, * 4th edition, McGraw-Hill, 1987, 407-410 * 5th edition, McGraw-Hill, 20001, p. 9.21/22 */ - Scalar muResult = 0; - const Scalar mu[numComponents] = { + LhsEval muResult = 0; + const LhsEval mu[numComponents] = { H2O::gasViscosity(T, H2O::vaporPressure(T)), N2::gasViscosity(T, p) }; - Scalar sumx = 0.0; + LhsEval sumx = 0.0; for (int compIdx = 0; compIdx < numComponents; ++compIdx) - sumx += fluidState.moleFraction(phaseIdx, compIdx); - sumx = std::max(1e-10, sumx); + sumx += FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, compIdx)); + sumx = LhsToolbox::max(1e-10, sumx); for (int i = 0; i < numComponents; ++i) { - Scalar divisor = 0; + LhsEval divisor = 0; for (int j = 0; j < numComponents; ++j) { - Scalar phiIJ = 1 + std::sqrt(mu[i]/mu[j]) * std::pow(molarMass(j)/molarMass(i), 1/4.0); + LhsEval phiIJ = 1 + LhsToolbox::sqrt(mu[i]/mu[j]) * std::pow(molarMass(j)/molarMass(i), 1/4.0); phiIJ *= phiIJ; phiIJ /= std::sqrt(8*(1 + molarMass(i)/molarMass(j))); - divisor += fluidState.moleFraction(phaseIdx, j)/sumx * phiIJ; + divisor += + FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, j)) + /sumx*phiIJ; } - muResult += fluidState.moleFraction(phaseIdx, i)/sumx * mu[i] / divisor; + muResult += + FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, i)) + /sumx*mu[i]/divisor; } return muResult; } } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); // liquid phase if (phaseIdx == liquidPhaseIdx) { @@ -390,21 +399,25 @@ public: return Opm::BinaryCoeff::H2O_N2::henry(T)/p; } + assert(phaseIdx == gasPhaseIdx); + // for the gas phase, assume an ideal gas when it comes to // fugacity (-> fugacity == partial pressure) return 1.0; } //! \copydoc BaseFluidSystem::diffusionCoefficient - template - static Scalar diffusionCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval diffusionCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); // liquid phase if (phaseIdx == liquidPhaseIdx) @@ -416,13 +429,15 @@ public: } //! \copydoc BaseFluidSystem::enthalpy - template - static Scalar enthalpy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); Valgrind::CheckDefined(T); Valgrind::CheckDefined(p); @@ -431,84 +446,86 @@ public: // TODO: correct way to deal with the solutes??? return H2O::liquidEnthalpy(T, p); } + // gas phase - else { - // assume ideal mixture: Molecules of one component don't - // "see" the molecules of the other component, which means - // that the total specific enthalpy is the sum of the - // "partial specific enthalpies" of the components. - Scalar hH2O = - fluidState.massFraction(gasPhaseIdx, H2OIdx) - * H2O::gasEnthalpy(T, p); - Scalar hN2 = - fluidState.massFraction(gasPhaseIdx, N2Idx) - * N2::gasEnthalpy(T, p); - return hH2O + hN2; - } + assert(phaseIdx == gasPhaseIdx); + + // assume ideal mixture: Molecules of one component don't + // "see" the molecules of the other component, which means + // that the total specific enthalpy is the sum of the + // "partial specific enthalpies" of the components. + const auto& XgH2O = FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, H2OIdx)); + const auto& XgN2 = FsToolbox::template toLhs(fluidState.massFraction(gasPhaseIdx, N2Idx)); + + LhsEval hH2O = XgH2O*H2O::gasEnthalpy(T, p); + LhsEval hN2 = XgN2*N2::gasEnthalpy(T, p); + return hH2O + hN2; } //! \copydoc BaseFluidSystem::thermalConductivity - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx) ; - Scalar pressure = fluidState.pressure(phaseIdx); + typedef Opm::MathToolbox FsToolbox; - if (phaseIdx == liquidPhaseIdx) { // liquid phase - return H2O::liquidThermalConductivity(temperature, pressure); - } - else { // gas phase - Scalar lambdaDryN2 = N2::gasThermalConductivity(temperature, pressure); - - if (useComplexRelations){ - Scalar xN2 = fluidState.moleFraction(phaseIdx, N2Idx); - Scalar xH2O = fluidState.moleFraction(phaseIdx, H2OIdx); - Scalar lambdaN2 = xN2 * lambdaDryN2; - - // Assuming Raoult's, Daltons law and ideal gas - // in order to obtain the partial density of water in the air phase - Scalar partialPressure = pressure * xH2O; - - Scalar lambdaH2O = - xH2O - * H2O::gasThermalConductivity(temperature, partialPressure); - return lambdaN2 + lambdaH2O; - } - else - return lambdaDryN2; // conductivity of Nitrogen [W / (m K ) ] + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + if (phaseIdx == liquidPhaseIdx) // liquid phase + return H2O::liquidThermalConductivity(T, p); + + // gas phase + assert(phaseIdx == gasPhaseIdx); + + if (useComplexRelations){ + // return the sum of the partial conductivity of Nitrogen and Steam + const auto& xH2O = FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, H2OIdx)); + const auto& xN2 = FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, N2Idx)); + + // Assuming Raoult's, Daltons law and ideal gas in order to obtain the + // partial pressures in the gas phase + const auto& lambdaN2 = N2::gasThermalConductivity(T, p*xN2); + const auto& lambdaH2O = H2O::gasThermalConductivity(T, p*xH2O); + + return lambdaN2 + lambdaH2O; } + else + // return the conductivity of dry Nitrogen + return N2::gasThermalConductivity(T, p); } //! \copydoc BaseFluidSystem::heatCapacity - template - static Scalar heatCapacity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval heatCapacity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { - if (phaseIdx == liquidPhaseIdx) { - return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx)); - } + typedef Opm::MathToolbox FsToolbox; + + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + const auto& xAlphaH2O = FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, H2OIdx)); + const auto& xAlphaN2 = FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, N2Idx)); + const auto& XAlphaH2O = FsToolbox::template toLhs(fluidState.massFraction(phaseIdx, H2OIdx)); + const auto& XAlphaN2 = FsToolbox::template toLhs(fluidState.massFraction(phaseIdx, N2Idx)); + + if (phaseIdx == liquidPhaseIdx) + return H2O::liquidHeatCapacity(T, p); + + assert(phaseIdx == gasPhaseIdx); // for the gas phase, assume ideal mixture, i.e. molecules of // one component don't "see" the molecules of the other // component - - Scalar c_pN2; - Scalar c_pH2O; + LhsEval c_pN2; + LhsEval c_pH2O; // let the water and nitrogen components do things their own way if (useComplexRelations) { - c_pN2 = N2::gasHeatCapacity(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx) - * fluidState.moleFraction(phaseIdx, N2Idx)); - - c_pH2O = H2O::gasHeatCapacity(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx) - * fluidState.moleFraction(phaseIdx, H2OIdx)); + c_pN2 = N2::gasHeatCapacity(T, p*xAlphaN2); + c_pH2O = H2O::gasHeatCapacity(T, p*xAlphaH2O); } else { // assume an ideal gas for both components. See: @@ -524,10 +541,9 @@ public: c_pH2O = c_pH2Omolar/molarMass(H2OIdx); } - // mangle both components together - return - c_pH2O*fluidState.massFraction(gasPhaseIdx, H2OIdx) - + c_pN2*fluidState.massFraction(gasPhaseIdx, N2Idx); + // mingle both components together. this assumes that there is no "cross + // interaction" between both flavors of molecules. + return XAlphaH2O*c_pH2O + XAlphaN2*c_pN2; } }; diff --git a/opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp b/opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp index 2c894bafc..2b969cb8b 100644 --- a/opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp +++ b/opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp @@ -245,19 +245,21 @@ public: } //! \copydoc BaseFluidSystem::density - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); - Scalar sumMoleFrac = 0; + LhsEval sumMoleFrac = 0; for (int compIdx = 0; compIdx < numComponents; ++compIdx) - sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx); + sumMoleFrac += FsToolbox::template toLhs(fluidState.moleFraction(phaseIdx, compIdx)); assert(phaseIdx == liquidPhaseIdx); @@ -267,47 +269,49 @@ public: else { // See: Ochs 2008 - Scalar rholH2O = H2O::liquidDensity(T, p); - Scalar clH2O = rholH2O/H2O::molarMass(); + const LhsEval& rholH2O = H2O::liquidDensity(T, p); + const LhsEval& clH2O = rholH2O/H2O::molarMass(); + + const auto& xlH2O = FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, H2OIdx)); + const auto& xlN2 = FsToolbox::template toLhs(fluidState.moleFraction(liquidPhaseIdx, N2Idx)); // this assumes each nitrogen molecule displaces exactly one // water molecule in the liquid - return - clH2O - * (H2O::molarMass()*fluidState.moleFraction(liquidPhaseIdx, H2OIdx) - + - N2::molarMass()*fluidState.moleFraction(liquidPhaseIdx, N2Idx)) - / sumMoleFrac; + return clH2O*(H2O::molarMass()*xlH2O + N2::molarMass()*xlN2)/sumMoleFrac; } } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(phaseIdx == liquidPhaseIdx); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); // assume pure water for the liquid phase return H2O::liquidViscosity(T, p); } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(phaseIdx == liquidPhaseIdx); assert(0 <= compIdx && compIdx < numComponents); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (compIdx == H2OIdx) return H2O::vaporPressure(T)/p; @@ -315,31 +319,35 @@ public: } //! \copydoc BaseFluidSystem::diffusionCoefficient - template - static Scalar diffusionCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval diffusionCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(phaseIdx == liquidPhaseIdx); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p); } //! \copydoc BaseFluidSystem::enthalpy - template - static Scalar enthalpy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert (phaseIdx == liquidPhaseIdx); - Scalar T = fluidState.temperature(phaseIdx); - Scalar p = fluidState.pressure(phaseIdx); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); Valgrind::CheckDefined(T); Valgrind::CheckDefined(p); @@ -348,32 +356,38 @@ public: } //! \copydoc BaseFluidSystem::thermalConductivity - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - const int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + const int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(phaseIdx == liquidPhaseIdx); if(useComplexRelations){ - Scalar temperature = fluidState.temperature(phaseIdx) ; - Scalar pressure = fluidState.pressure(phaseIdx); - return H2O::liquidThermalConductivity(temperature, pressure); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + return H2O::liquidThermalConductivity(T, p); } else return 0.578078; // conductivity of water[W / (m K ) ] IAPWS evaluated at p=.1 MPa, T=8C } //! \copydoc BaseFluidSystem::heatCapacity - template - static Scalar heatCapacity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval heatCapacity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert (phaseIdx == liquidPhaseIdx); - return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx), - fluidState.pressure(phaseIdx)); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + + return H2O::liquidHeatCapacity(T, p); } }; diff --git a/opm/material/fluidsystems/LiquidPhase.hpp b/opm/material/fluidsystems/LiquidPhase.hpp index 9fe1dcf27..5dcac1134 100644 --- a/opm/material/fluidsystems/LiquidPhase.hpp +++ b/opm/material/fluidsystems/LiquidPhase.hpp @@ -74,35 +74,43 @@ public: { return Component::triplePressure(); } //! \copydoc GasPhase::vaporPressure - static Scalar vaporPressure(Scalar temperature) + template + static Evaluation vaporPressure(const Evaluation& temperature) { return Component::vaporPressure(temperature); } //! \copydoc GasPhase::density - static Scalar density(Scalar temperature, Scalar pressure) + template + static Evaluation density(const Evaluation& temperature, const Evaluation& pressure) { return Component::liquidDensity(temperature, pressure); } //! \copydoc GasPhase::pressure - static Scalar pressure(Scalar temperature, Scalar density) + template + static Evaluation pressure(const Evaluation& temperature, const Evaluation& density) { return Component::liquidPressure(temperature, density); } //! \copydoc GasPhase::enthalpy - static const Scalar enthalpy(Scalar temperature, Scalar pressure) + template + static const Evaluation enthalpy(const Evaluation& temperature, const Evaluation& pressure) { return Component::liquidEnthalpy(temperature, pressure); } //! \copydoc GasPhase::internalEnergy - static const Scalar internalEnergy(Scalar temperature, Scalar pressure) + template + static const Evaluation internalEnergy(const Evaluation& temperature, const Evaluation& pressure) { return Component::liquidInternalEnergy(temperature, pressure); } //! \copydoc GasPhase::viscosity - static Scalar viscosity(Scalar temperature, Scalar pressure) + template + static Evaluation viscosity(const Evaluation& temperature, const Evaluation& pressure) { return Component::liquidViscosity(temperature, pressure); } //! \copydoc GasPhase::thermalConductivity - static Scalar thermalConductivity(Scalar temperature, Scalar pressure) + template + static Evaluation thermalConductivity(const Evaluation& temperature, const Evaluation& pressure) { return Component::liquidThermalConductivity(temperature, pressure); } //! \copydoc GasPhase::heatCapacity - static Scalar heatCapacity(Scalar temperature, Scalar pressure) + template + static Evaluation heatCapacity(const Evaluation& temperature, const Evaluation& pressure) { return Component::liquidHeatCapacity(temperature, pressure); } }; } // namespace Opm diff --git a/opm/material/fluidsystems/SinglePhaseFluidSystem.hpp b/opm/material/fluidsystems/SinglePhaseFluidSystem.hpp index 2bb806176..77245e02a 100644 --- a/opm/material/fluidsystems/SinglePhaseFluidSystem.hpp +++ b/opm/material/fluidsystems/SinglePhaseFluidSystem.hpp @@ -126,7 +126,7 @@ public: } //! \copydoc BaseFluidSystem::molarMass - static const Scalar molarMass(int compIdx) + static Scalar molarMass(int compIdx) { //assert(0 <= compIdx && compIdx < numComponents); @@ -138,7 +138,7 @@ public: * * \param compIdx The index of the component to consider */ - static const Scalar criticalTemperature(int compIdx) + static Scalar criticalTemperature(int compIdx) { //assert(0 <= compIdx && compIdx < numComponents); @@ -150,7 +150,7 @@ public: * * \param compIdx The index of the component to consider */ - static const Scalar criticalPressure(int compIdx) + static Scalar criticalPressure(int compIdx) { //assert(0 <= compIdx && compIdx < numComponents); @@ -162,7 +162,7 @@ public: * * \param compIdx The index of the component to consider */ - static const Scalar acentricFactor(int compIdx) + static Scalar acentricFactor(int compIdx) { //assert(0 <= compIdx && compIdx < numComponents); @@ -178,37 +178,41 @@ public: { } //! \copydoc BaseFluidSystem::density - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); - return Fluid::density(temperature, pressure); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + return Fluid::density(T, p); } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); - return Fluid::viscosity(temperature, pressure); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + return Fluid::viscosity(T, p); } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); @@ -223,42 +227,48 @@ public: } //! \copydoc BaseFluidSystem::enthalpy - template - static Scalar enthalpy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); - return Fluid::enthalpy(temperature, pressure); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + return Fluid::enthalpy(T, p); } //! \copydoc BaseFluidSystem::thermalConductivity - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); - return Fluid::thermalConductivity(temperature, pressure); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + return Fluid::thermalConductivity(T, p); } //! \copydoc BaseFluidSystem::heatCapacity - template - static Scalar heatCapacity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval heatCapacity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef Opm::MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); - return Fluid::heatCapacity(temperature, pressure); + const auto& T = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& p = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); + return Fluid::heatCapacity(T, p); } }; diff --git a/opm/material/fluidsystems/Spe5FluidSystem.hpp b/opm/material/fluidsystems/Spe5FluidSystem.hpp index 5c9bbea2f..ed3f8b72d 100644 --- a/opm/material/fluidsystems/Spe5FluidSystem.hpp +++ b/opm/material/fluidsystems/Spe5FluidSystem.hpp @@ -352,23 +352,27 @@ public: } //! \copydoc BaseFluidSystem::density - template + template static Scalar density(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { assert(0 <= phaseIdx && phaseIdx < numPhases); + static_assert(std::is_same::value, + "The SPE-5 fluid system is currently only implemented for the scalar case."); return fluidState.averageMolarMass(phaseIdx)/paramCache.molarVolume(phaseIdx); } //! \copydoc BaseFluidSystem::viscosity - template + template static Scalar viscosity(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx) { assert(0 <= phaseIdx && phaseIdx <= numPhases); + static_assert(std::is_same::value, + "The SPE-5 fluid system is currently only implemented for the scalar case."); if (phaseIdx == gasPhaseIdx) { // given by SPE-5 in table on page 64. we use a constant @@ -387,7 +391,7 @@ public: } //! \copydoc BaseFluidSystem::fugacityCoefficient - template + template static Scalar fugacityCoefficient(const FluidState &fluidState, const ParameterCache ¶mCache, int phaseIdx, @@ -395,6 +399,8 @@ public: { assert(0 <= phaseIdx && phaseIdx <= numPhases); assert(0 <= compIdx && compIdx <= numComponents); + static_assert(std::is_same::value, + "The SPE-5 fluid system is currently only implemented for the scalar case."); if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx) return PengRobinsonMixture::computeFugacityCoefficient(fluidState, diff --git a/opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp b/opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp index 43d201366..b0e0d0c77 100644 --- a/opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp +++ b/opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp @@ -52,7 +52,7 @@ namespace FluidSystems { */ template class TwoPhaseImmiscible -: public BaseFluidSystem > + : public BaseFluidSystem > { // do not try to instanciate this class, it has only static members! TwoPhaseImmiscible() @@ -216,42 +216,48 @@ public: } //! \copydoc BaseFluidSystem::density - template - static Scalar density(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval density(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const auto& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == wettingPhaseIdx) return WettingPhase::density(temperature, pressure); return NonwettingPhase::density(temperature, pressure); } //! \copydoc BaseFluidSystem::viscosity - template - static Scalar viscosity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval viscosity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const auto& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == wettingPhaseIdx) return WettingPhase::viscosity(temperature, pressure); return NonwettingPhase::viscosity(temperature, pressure); } //! \copydoc BaseFluidSystem::fugacityCoefficient - template - static Scalar fugacityCoefficient(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx, - int compIdx) + template + static LhsEval fugacityCoefficient(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx, + int compIdx) { + typedef MathToolbox LhsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); assert(0 <= compIdx && compIdx < numComponents); @@ -260,50 +266,56 @@ public: // the component in the fluid. Probably that's not worth // the effort, since the fugacity coefficient of the other // component is infinite anyway... - return 1.0; - return std::numeric_limits::infinity(); + return LhsToolbox::createConstant(1.0); + return LhsToolbox::createConstant(std::numeric_limits::infinity()); } //! \copydoc BaseFluidSystem::enthalpy - template - static Scalar enthalpy(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval enthalpy(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const auto& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == wettingPhaseIdx) return WettingPhase::enthalpy(temperature, pressure); return NonwettingPhase::enthalpy(temperature, pressure); } //! \copydoc BaseFluidSystem::thermalConductivity - template - static Scalar thermalConductivity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval thermalConductivity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const auto& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == wettingPhaseIdx) return WettingPhase::thermalConductivity(temperature, pressure); return NonwettingPhase::thermalConductivity(temperature, pressure); } //! \copydoc BaseFluidSystem::heatCapacity - template - static Scalar heatCapacity(const FluidState &fluidState, - const ParameterCache ¶mCache, - int phaseIdx) + template + static LhsEval heatCapacity(const FluidState &fluidState, + const ParameterCache ¶mCache, + int phaseIdx) { + typedef MathToolbox FsToolbox; + assert(0 <= phaseIdx && phaseIdx < numPhases); - Scalar temperature = fluidState.temperature(phaseIdx); - Scalar pressure = fluidState.pressure(phaseIdx); + const auto& temperature = FsToolbox::template toLhs(fluidState.temperature(phaseIdx)); + const auto& pressure = FsToolbox::template toLhs(fluidState.pressure(phaseIdx)); if (phaseIdx == wettingPhaseIdx) return WettingPhase::heatCapacity(temperature, pressure); return NonwettingPhase::heatCapacity(temperature, pressure); diff --git a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp index 47c74148f..097578c16 100644 --- a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp @@ -40,10 +40,16 @@ namespace Opm { * \brief This class represents the Pressure-Volume-Temperature relations of the oil phase * without dissolved gas and constant compressibility/"viscosibility". */ -template -class ConstantCompressibilityOilPvt : public OilPvtInterface +template +class ConstantCompressibilityOilPvt : public OilPvtInterfaceTemplateWrapper > { - typedef FluidSystems::BlackOil BlackOilFluidSystem; + friend class OilPvtInterfaceTemplateWrapper >; + + typedef FluidSystems::BlackOil BlackOilFluidSystem; typedef Opm::Tabulated1DFunction TabulatedOneDFunction; typedef std::vector > SamplingPoints; @@ -133,22 +139,24 @@ public: void initEnd() { } +private: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ - Scalar viscosity(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval viscosity_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { // Eclipse calculates the viscosity in a weird way: it // calcultes the product of B_w and mu_w and then divides the // result by B_w... Scalar BoMuoRef = oilViscosity_[regionIdx]*oilReferenceFormationVolumeFactor_[regionIdx]; - Scalar Bo = formationVolumeFactor(regionIdx, temperature, pressure, XoG); + const LhsEval& Bo = formationVolumeFactor_(regionIdx, temperature, pressure, XoG); Scalar pRef = oilReferencePressure_[regionIdx]; - Scalar Y = + const LhsEval& Y = (oilCompressibility_[regionIdx] - oilViscosibility_[regionIdx]) * (pressure - pRef); return BoMuoRef/((1 + Y*(1 + Y/2))*Bo); @@ -157,12 +165,13 @@ public: /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ - Scalar density(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval density_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { - Scalar Bo = formationVolumeFactor(regionIdx, temperature, pressure, XoG); + const LhsEval& Bo = formationVolumeFactor_(regionIdx, temperature, pressure, XoG); Scalar rhooRef = BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx); return rhooRef/Bo; } @@ -170,14 +179,15 @@ public: /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ - Scalar formationVolumeFactor(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval formationVolumeFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { // cf. ECLiPSE 2011 technical description, p. 116 Scalar pRef = oilReferencePressure_[regionIdx]; - Scalar X = oilCompressibility_[regionIdx]*(pressure - pRef); + const LhsEval& X = oilCompressibility_[regionIdx]*(pressure - pRef); Scalar BoRef = oilReferenceFormationVolumeFactor_[regionIdx]; return BoRef/(1 + X*(1 + X/2)); @@ -187,15 +197,16 @@ public: * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ - Scalar fugacityCoefficient(int regionIdx, - Scalar temperature, - Scalar pressure, - int compIdx) const OPM_FINAL + template + LhsEval fugacityCoefficient_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int compIdx) const { // set the oil component fugacity coefficient in oil phase // arbitrarily. we use some pseudo-realistic value for the vapor // pressure to ease physical interpretation of the results - Scalar phi_oO = 20e3/pressure; + const LhsEval& phi_oO = 20e3/pressure; if (compIdx == BlackOilFluidSystem::oilCompIdx) return phi_oO; @@ -213,10 +224,15 @@ public: /*! * \brief Returns the gas dissolution factor \f$R_s\f$ [m^3/m^3] of the oil phase. */ - Scalar gasDissolutionFactor(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; /* this is dead oil! */ } + template + LhsEval gasDissolutionFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + + return Toolbox::createConstant(0.0); /* this is dead oil! */ + } /*! * \brief Returns the saturation pressure of the oil phase [Pa] @@ -224,20 +240,35 @@ public: * * \param XoG The mass fraction of the gas component in the oil phase [-] */ - Scalar oilSaturationPressure(int regionIdx, - Scalar temperature, - Scalar XoG) const OPM_FINAL - { return 0.0; /* this is dead oil, so there isn't any meaningful saturation pressure! */ } + template + LhsEval oilSaturationPressure_(int regionIdx, + const LhsEval& temperature, + const LhsEval& XoG) const + { + typedef Opm::MathToolbox Toolbox; - Scalar saturatedOilGasMassFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; /* this is dead oil! */ } + return Toolbox::createConstant(0.0); /* this is dead oil, so there isn't any meaningful saturation pressure! */ + } - Scalar saturatedOilGasMoleFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; /* this is dead oil! */ } + template + LhsEval saturatedOilGasMassFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + + return Toolbox::createConstant(0.0); /* this is dead oil! */ + } + + template + LhsEval saturatedOilGasMoleFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + + return Toolbox::createConstant(0.0); /* this is dead oil! */ + } private: std::vector oilReferencePressure_; diff --git a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp index b6803ea3b..9e5554dcf 100644 --- a/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp @@ -26,7 +26,6 @@ #include "WaterPvtInterface.hpp" #include - #if HAVE_OPM_PARSER #include #endif @@ -39,10 +38,13 @@ namespace Opm { * \brief This class represents the Pressure-Volume-Temperature relations of the gas phase * without vaporized oil. */ -template -class ConstantCompressibilityWaterPvt : public WaterPvtInterface +template +class ConstantCompressibilityWaterPvt + : public WaterPvtInterfaceTemplateWrapper > { - typedef FluidSystems::BlackOil BlackOilFluidSystem; + friend class WaterPvtInterfaceTemplateWrapper >; + + typedef FluidSystems::BlackOil BlackOilFluidSystem; typedef Opm::Tabulated1DFunction TabulatedOneDFunction; typedef std::vector > SamplingPoints; @@ -132,21 +134,23 @@ public: void initEnd() { } +private: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ - Scalar viscosity(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL + template + LhsEval viscosity_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const { // Eclipse calculates the viscosity in a weird way: it // calcultes the product of B_w and mu_w and then divides the // result by B_w... Scalar BwMuwRef = waterViscosity_[regionIdx]*waterReferenceFormationVolumeFactor_[regionIdx]; - Scalar Bw = formationVolumeFactor(regionIdx, temperature, pressure); + const LhsEval& Bw = formationVolumeFactor_(regionIdx, temperature, pressure); Scalar pRef = waterReferencePressure_[regionIdx]; - Scalar Y = + const LhsEval& Y = (waterCompressibility_[regionIdx] - waterViscosibility_[regionIdx]) * (pressure - pRef); return BwMuwRef/((1 + Y*(1 + Y/2))*Bw); @@ -155,11 +159,12 @@ public: /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ - Scalar density(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL + template + LhsEval density_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const { - Scalar Bw = formationVolumeFactor(regionIdx, temperature, pressure); + const LhsEval& Bw = formationVolumeFactor_(regionIdx, temperature, pressure); Scalar rhowRef = BlackOilFluidSystem::referenceDensity(waterPhaseIdx, regionIdx); return rhowRef/Bw; } @@ -167,13 +172,14 @@ public: /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ - Scalar formationVolumeFactor(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL + template + LhsEval formationVolumeFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const { // cf. ECLiPSE 2011 technical description, p. 116 Scalar pRef = waterReferencePressure_[regionIdx]; - Scalar X = waterCompressibility_[regionIdx]*(pressure - pRef); + const LhsEval& X = waterCompressibility_[regionIdx]*(pressure - pRef); Scalar BwRef = waterReferenceFormationVolumeFactor_[regionIdx]; @@ -185,10 +191,11 @@ public: * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ - Scalar fugacityCoefficient(int regionIdx, - Scalar temperature, - Scalar pressure, - int compIdx) const OPM_FINAL + template + LhsEval fugacityCoefficient_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int compIdx) const { // set the affinity of the gas and oil components to the water phase to be 10 // orders of magnitute smaller than that of the water component. for this we use diff --git a/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp b/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp index 160f33c29..95d44c9a8 100644 --- a/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp @@ -40,10 +40,16 @@ namespace Opm { * \brief This class represents the Pressure-Volume-Temperature relations of the oil phase * without dissolved gas. */ -template -class DeadOilPvt : public OilPvtInterface +template +class DeadOilPvt : public OilPvtInterfaceTemplateWrapper > { - typedef FluidSystems::BlackOil BlackOilFluidSystem; + friend class OilPvtInterfaceTemplateWrapper >; + + typedef FluidSystems::BlackOil BlackOilFluidSystem; typedef Opm::Tabulated1DFunction TabulatedOneDFunction; typedef std::vector > SamplingPoints; @@ -96,8 +102,8 @@ public: { assert(pvdoTable.numRows() > 1); - const std::vector& BColumn(pvdoTable.getFormationFactorColumn()); - std::vector invBColumn(pvdoTable.getFormationFactorColumn()); + const auto& BColumn(pvdoTable.getFormationFactorColumn()); + std::vector invBColumn(BColumn.size()); for (unsigned i = 0; i < invBColumn.size(); ++i) invBColumn[i] = 1/BColumn[i]; @@ -141,16 +147,18 @@ public: } } +private: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ - Scalar viscosity(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval viscosity_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { - Scalar invBo = inverseOilB_[regionIdx].eval(pressure, /*extrapolate=*/true); - Scalar invMuoBo = inverseOilBMu_[regionIdx].eval(pressure, /*extrapolate=*/true); + const LhsEval& invBo = inverseOilB_[regionIdx].eval(pressure, /*extrapolate=*/true); + const LhsEval& invMuoBo = inverseOilBMu_[regionIdx].eval(pressure, /*extrapolate=*/true); return invBo/invMuoBo; } @@ -158,39 +166,42 @@ public: /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ - Scalar density(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval density_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { Scalar rhooRef = BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx); - Scalar Bo = formationVolumeFactor(regionIdx, temperature, pressure, XoG); + const LhsEval& Bo = formationVolumeFactor_(regionIdx, temperature, pressure, XoG); return rhooRef/Bo; } /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ - Scalar formationVolumeFactor(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval formationVolumeFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { return 1.0 / inverseOilB_[regionIdx].eval(pressure, /*extrapolate=*/true); } /*! * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ - Scalar fugacityCoefficient(int regionIdx, - Scalar temperature, - Scalar pressure, - int compIdx) const OPM_FINAL + template + LhsEval fugacityCoefficient_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int compIdx) const { // set the oil component fugacity coefficient in oil phase // arbitrarily. we use some pseudo-realistic value for the vapor // pressure to ease physical interpretation of the results - Scalar phi_oO = 20e3/pressure; + const LhsEval& phi_oO = 20e3/pressure; if (compIdx == BlackOilFluidSystem::oilCompIdx) return phi_oO; @@ -208,10 +219,15 @@ public: /*! * \brief Returns the gas dissolution factor \f$R_s\f$ [m^3/m^3] of the oil phase. */ - Scalar gasDissolutionFactor(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; /* this is dead oil! */ } + template + LhsEval gasDissolutionFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + + return Toolbox::createConstant(0.0); /* this is dead oil! */ + } /*! * \brief Returns the saturation pressure of the oil phase [Pa] @@ -219,20 +235,35 @@ public: * * \param XoG The mass fraction of the gas component in the oil phase [-] */ - Scalar oilSaturationPressure(int regionIdx, - Scalar temperature, - Scalar XoG) const OPM_FINAL - { return 0.0; /* this is dead oil, so there isn't any meaningful saturation pressure! */ } + template + LhsEval oilSaturationPressure_(int regionIdx, + const LhsEval& temperature, + const LhsEval& XoG) const + { + typedef Opm::MathToolbox Toolbox; - Scalar saturatedOilGasMassFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; /* this is dead oil! */ } + return Toolbox::createConstant(0.0); /* this is dead oil, so there isn't any meaningful saturation pressure! */ + } - Scalar saturatedOilGasMoleFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; /* this is dead oil! */ } + template + LhsEval saturatedOilGasMassFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + + return Toolbox::createConstant(0.0); /* this is dead oil! */ + } + + template + LhsEval saturatedOilGasMoleFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + + return Toolbox::createConstant(0.0); /* this is dead oil! */ + } private: std::vector inverseOilB_; diff --git a/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp b/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp index 000018fa6..c91e9cca8 100644 --- a/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp @@ -40,10 +40,13 @@ namespace Opm { * \brief This class represents the Pressure-Volume-Temperature relations of the gas phase * without vaporized oil. */ -template -class DryGasPvt : public GasPvtInterface +template +class DryGasPvt + : public GasPvtInterfaceTemplateWrapper > { - typedef FluidSystems::BlackOil BlackOilFluidSystem; + friend class GasPvtInterfaceTemplateWrapper >; + + typedef FluidSystems::BlackOil BlackOilFluidSystem; typedef Opm::Tabulated1DFunction TabulatedOneDFunction; typedef std::vector > SamplingPoints; @@ -133,16 +136,18 @@ public: } } +private: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ - Scalar viscosity(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XgO) const OPM_FINAL + template + LhsEval viscosity_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO) const { - Scalar invBg = inverseGasB_[regionIdx].eval(pressure, /*extrapolate=*/true); - Scalar invMugBg = inverseGasBMu_[regionIdx].eval(pressure, /*extrapolate=*/true); + const LhsEval& invBg = inverseGasB_[regionIdx].eval(pressure, /*extrapolate=*/true); + const LhsEval& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, /*extrapolate=*/true); return invBg/invMugBg; } @@ -150,38 +155,43 @@ public: /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ - Scalar density(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XgO) const OPM_FINAL + template + LhsEval density_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO) const { // gas formation volume factor at reservoir pressure - Scalar Bg = formationVolumeFactor(regionIdx, temperature, pressure, XgO); + const LhsEval& Bg = formationVolumeFactor_(regionIdx, temperature, pressure, XgO); return BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx)/Bg; } /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ - Scalar formationVolumeFactor(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XgO) const OPM_FINAL + template + LhsEval formationVolumeFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO) const { return 1.0/inverseGasB_[regionIdx].eval(pressure, /*extrapolate=*/true); } /*! * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ - Scalar fugacityCoefficient(int regionIdx, - Scalar temperature, - Scalar pressure, - int compIdx) const OPM_FINAL + template + LhsEval fugacityCoefficient_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int compIdx) const { + typedef Opm::MathToolbox Toolbox; + // make the gas component more affine to the gas phase than the other components if (compIdx == BlackOilFluidSystem::gasCompIdx) - return 1; - return 1e6; + return Toolbox::createConstant(1.0); + return Toolbox::createConstant(1e6); } /*! @@ -190,28 +200,44 @@ public: * * \param XgO The mass fraction of the oil component in the gas phase [-] */ - Scalar gasSaturationPressure(int regionIdx, - Scalar temperature, - Scalar XgO) const OPM_FINAL - { return 0.0; } // this is dry gas! + template + LhsEval gasSaturationPressure_(int regionIdx, + const LhsEval& temperature, + const LhsEval& XgO) const + { + typedef Opm::MathToolbox Toolbox; + return Toolbox::createConstant(0.0); // this is dry gas! + } /*! * \brief Returns the gas dissolution factor \f$R_s\f$ [m^3/m^3] of the oil phase. */ - Scalar oilVaporizationFactor(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; } // this is dry gas! + template + LhsEval oilVaporizationFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + return Toolbox::createConstant(0.0); // this is dry gas! + } - Scalar saturatedGasOilMassFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; } // this is dry gas! + template + LhsEval saturatedGasOilMassFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + return Toolbox::createConstant(0.0); // this is dry gas! + } - Scalar saturatedGasOilMoleFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return 0.0; } // this is dry gas! + template + LhsEval saturatedGasOilMoleFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { + typedef Opm::MathToolbox Toolbox; + return Toolbox::createConstant(0.0); // this is dry gas! + } private: std::vector inverseGasB_; diff --git a/opm/material/fluidsystems/blackoilpvt/GasPvtInterface.hpp b/opm/material/fluidsystems/blackoilpvt/GasPvtInterface.hpp index fc016a3dd..5463cad5b 100644 --- a/opm/material/fluidsystems/blackoilpvt/GasPvtInterface.hpp +++ b/opm/material/fluidsystems/blackoilpvt/GasPvtInterface.hpp @@ -23,8 +23,9 @@ #ifndef OPM_GAS_PVT_INTERFACE_HPP #define OPM_GAS_PVT_INTERFACE_HPP -namespace Opm { +#include +namespace Opm { /*! * \brief This class represents the Pressure-Volume-Temperature relations of the gas * phase in the black-oil model. @@ -37,38 +38,54 @@ namespace Opm { * Note that, since the application for this class is the black oil fluid system, the API * exposed by this class is pretty specific to the black oil model. */ -template +template class GasPvtInterface { public: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ + virtual Evaluation viscosity(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XgO) const = 0; virtual Scalar viscosity(int regionIdx, Scalar temperature, Scalar pressure, - Scalar XoG) const = 0; + Scalar XgO) const = 0; /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ + virtual Evaluation formationVolumeFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XgO) const = 0; virtual Scalar formationVolumeFactor(int regionIdx, Scalar temperature, Scalar pressure, - Scalar XoG) const = 0; + Scalar XgO) const = 0; /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ + virtual Evaluation density(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XgO) const = 0; virtual Scalar density(int regionIdx, Scalar temperature, Scalar pressure, - Scalar XoG) const = 0; + Scalar XgO) const = 0; /*! * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ + virtual Evaluation fugacityCoefficient(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + int compIdx) const = 0; virtual Scalar fugacityCoefficient(int regionIdx, Scalar temperature, Scalar pressure, @@ -77,6 +94,9 @@ public: /*! * \brief Returns the oil vaporization factor \f$R_v\f$ [m^3/m^3] of oil saturated gas. */ + virtual Evaluation oilVaporizationFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; virtual Scalar oilVaporizationFactor(int regionIdx, Scalar temperature, Scalar pressure) const = 0; @@ -87,6 +107,9 @@ public: * * \param XgO The mass fraction of the oil component in the gas phase [-] */ + virtual Evaluation gasSaturationPressure(int regionIdx, + const Evaluation& temperature, + const Evaluation& XgO) const = 0; virtual Scalar gasSaturationPressure(int regionIdx, Scalar temperature, Scalar XgO) const = 0; @@ -95,6 +118,9 @@ public: * \brief Returns the gas mass fraction of oil-saturated gas at a given temperatire * and pressure [-]. */ + virtual Evaluation saturatedGasOilMassFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; virtual Scalar saturatedGasOilMassFraction(int regionIdx, Scalar temperature, Scalar pressure) const = 0; @@ -103,11 +129,182 @@ public: * \brief Returns the gas mole fraction of oil-saturated gas at a given temperatire * and pressure [-]. */ + virtual Evaluation saturatedGasOilMoleFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; virtual Scalar saturatedGasOilMoleFraction(int regionIdx, Scalar temperature, Scalar pressure) const = 0; }; +template +class GasPvtInterface +{ +public: + virtual Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const = 0; + virtual Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const = 0; + virtual Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const = 0; + virtual Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const = 0; + virtual Scalar oilVaporizationFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; + virtual Scalar gasSaturationPressure(int regionIdx, + Scalar temperature, + Scalar XgO) const = 0; + virtual Scalar saturatedGasOilMassFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; + virtual Scalar saturatedGasOilMoleFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; +}; + +template +class GasPvtInterfaceTemplateWrapper : public GasPvtInterface +{ + const Implementation& asImp_() const + { return *static_cast(this); } + +public: + virtual Evaluation viscosity(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XgO) const OPM_FINAL + { return asImp_().viscosity_(regionIdx, temperature, pressure, XgO); }; + virtual Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const OPM_FINAL + { return asImp_().viscosity_(regionIdx, temperature, pressure, XgO); }; + + virtual Evaluation formationVolumeFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XgO) const OPM_FINAL + { return asImp_().formationVolumeFactor_(regionIdx, temperature, pressure, XgO); }; + virtual Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const OPM_FINAL + { return asImp_().formationVolumeFactor_(regionIdx, temperature, pressure, XgO); }; + + virtual Evaluation density(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XgO) const OPM_FINAL + { return asImp_().density_(regionIdx, temperature, pressure, XgO); }; + virtual Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const OPM_FINAL + { return asImp_().density_(regionIdx, temperature, pressure, XgO); }; + + virtual Evaluation fugacityCoefficient(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + int compIdx) const OPM_FINAL + { return asImp_().fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); }; + virtual Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const OPM_FINAL + { return asImp_().fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); }; + + virtual Evaluation oilVaporizationFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().oilVaporizationFactor_(regionIdx, temperature, pressure); }; + virtual Scalar oilVaporizationFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().oilVaporizationFactor_(regionIdx, temperature, pressure); }; + + virtual Evaluation gasSaturationPressure(int regionIdx, + const Evaluation& temperature, + const Evaluation& XgO) const OPM_FINAL + { return asImp_().gasSaturationPressure_(regionIdx, temperature, XgO); }; + virtual Scalar gasSaturationPressure(int regionIdx, + Scalar temperature, + Scalar XgO) const OPM_FINAL + { return asImp_().gasSaturationPressure_(regionIdx, temperature, XgO); }; + + virtual Evaluation saturatedGasOilMassFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().saturatedGasOilMassFraction_(regionIdx, temperature, pressure); }; + virtual Scalar saturatedGasOilMassFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().saturatedGasOilMassFraction_(regionIdx, temperature, pressure); }; + + virtual Evaluation saturatedGasOilMoleFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().saturatedGasOilMoleFraction_(regionIdx, temperature, pressure); }; + virtual Scalar saturatedGasOilMoleFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().saturatedGasOilMoleFraction_(regionIdx, temperature, pressure); }; +}; + +template +class GasPvtInterfaceTemplateWrapper + : public GasPvtInterface +{ + const Implementation& asImp_() const + { return *static_cast(this); } + +public: + virtual Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const OPM_FINAL + { return asImp_().viscosity_(regionIdx, temperature, pressure, XgO); }; + virtual Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const OPM_FINAL + { return asImp_().formationVolumeFactor_(regionIdx, temperature, pressure, XgO); }; + virtual Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XgO) const OPM_FINAL + { return asImp_().density_(regionIdx, temperature, pressure, XgO); }; + virtual Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const OPM_FINAL + { return asImp_().fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); }; + virtual Scalar oilVaporizationFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().oilVaporizationFactor_(regionIdx, temperature, pressure); }; + virtual Scalar gasSaturationPressure(int regionIdx, + Scalar temperature, + Scalar XgO) const OPM_FINAL + { return asImp_().gasSaturationPressure_(regionIdx, temperature, XgO); }; + virtual Scalar saturatedGasOilMassFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().saturatedGasOilMassFraction_(regionIdx, temperature, pressure); }; + virtual Scalar saturatedGasOilMoleFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().saturatedGasOilMoleFraction_(regionIdx, temperature, pressure); }; +}; + } // namespace Opm #endif diff --git a/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp b/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp index 8036693df..be076fa8b 100644 --- a/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp @@ -40,10 +40,12 @@ namespace Opm { * \brief This class represents the Pressure-Volume-Temperature relations of the oil phas * with dissolved gas. */ -template -class LiveOilPvt : public OilPvtInterface +template +class LiveOilPvt : public OilPvtInterfaceTemplateWrapper > { - typedef FluidSystems::BlackOil BlackOilFluidSystem; + friend class OilPvtInterfaceTemplateWrapper >; + + typedef FluidSystems::BlackOil BlackOilFluidSystem; typedef Opm::UniformXTabulated2DFunction TabulatedTwoDFunction; typedef Opm::Tabulated1DFunction TabulatedOneDFunction; @@ -61,11 +63,11 @@ class LiveOilPvt : public OilPvtInterface public: void setNumRegions(int numRegions) { - if (static_cast(inverseOilB_.size()) < numRegions) { - inverseOilB_.resize(numRegions); - inverseOilBMu_.resize(numRegions); - oilMu_.resize(numRegions); - gasDissolutionFactor_.resize(numRegions); + if (static_cast(inverseOilBTable_.size()) < numRegions) { + inverseOilBTable_.resize(numRegions); + inverseOilBMuTable_.resize(numRegions); + oilMuTable_.resize(numRegions); + gasDissolutionFactorTable_.resize(numRegions); saturationPressureSpline_.resize(numRegions); } } @@ -76,7 +78,7 @@ public: * \param samplePoints A container of (x,y) values. */ void setSaturatedOilGasDissolutionFactor(int regionIdx, const SamplingPoints &samplePoints) - { gasDissolutionFactor_[regionIdx].setContainerOfTuples(samplePoints); } + { gasDissolutionFactorTable_[regionIdx].setContainerOfTuples(samplePoints); } /*! * \brief Initialize the function for the oil formation volume factor @@ -89,14 +91,14 @@ public: */ void setSaturatedOilFormationVolumeFactor(int regionIdx, const SamplingPoints &samplePoints) { - auto& invOilB = inverseOilB_[regionIdx]; + auto& invOilB = inverseOilBTable_[regionIdx]; - auto &Rs = gasDissolutionFactor_[regionIdx]; + auto &Rs = gasDissolutionFactorTable_[regionIdx]; Scalar T = BlackOilFluidSystem::surfaceTemperature; Scalar RsMin = 0.0; - Scalar RsMax = Rs.eval(gasDissolutionFactor_[regionIdx].xMax(), /*extrapolate=*/true); + Scalar RsMax = Rs.eval(gasDissolutionFactorTable_[regionIdx].xMax(), /*extrapolate=*/true); Scalar poMin = samplePoints.front().first; Scalar poMax = samplePoints.back().first; @@ -123,7 +125,7 @@ public: for (size_t pIdx = 0; pIdx < nP; ++pIdx) { Scalar po = poMin + (poMax - poMin)*pIdx/nP; - Scalar poSat = oilSaturationPressure(regionIdx, T, XoG); + Scalar poSat = oilSaturationPressure_(regionIdx, T, XoG); Scalar BoSat = oilFormationVolumeFactorSpline.eval(poSat, /*extrapolate=*/true); Scalar drhoo_dp = (1.1200 - 1.1189)/((5000 - 4000)*6894.76); Scalar rhoo = BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx)/BoSat*(1 + drhoo_dp*(po - poSat)); @@ -148,7 +150,7 @@ public: * and not the other way around. */ void setInverseOilFormationVolumeFactor(int regionIdx, const TabulatedTwoDFunction& invBo) - { inverseOilB_[regionIdx] = invBo; } + { inverseOilBTable_[regionIdx] = invBo; } /*! * \brief Initialize the viscosity of the oil phase. @@ -156,7 +158,7 @@ public: * This is a function of \f$(R_s, p_o)\f$... */ void setOilViscosity(int regionIdx, const TabulatedTwoDFunction& muo) - { oilMu_[regionIdx] = muo; } + { oilMuTable_[regionIdx] = muo; } /*! * \brief Initialize the phase viscosity for gas saturated oil @@ -167,10 +169,10 @@ public: */ void setSaturatedOilViscosity(int regionIdx, const SamplingPoints &samplePoints ) { - auto& gasDissolutionFactor = gasDissolutionFactor_[regionIdx]; + auto& gasDissolutionFactor = gasDissolutionFactorTable_[regionIdx]; Scalar RsMin = 0.0; - Scalar RsMax = gasDissolutionFactor.eval(gasDissolutionFactor_[regionIdx].xMax(), /*extrapolate=*/true); + Scalar RsMax = gasDissolutionFactor.eval(gasDissolutionFactorTable_[regionIdx].xMax(), /*extrapolate=*/true); Scalar poMin = samplePoints.front().first; Scalar poMax = samplePoints.back().first; @@ -186,13 +188,13 @@ public: for (size_t RsIdx = 0; RsIdx < nRs; ++RsIdx) { Scalar Rs = RsMin + (RsMax - RsMin)*RsIdx/nRs; - oilMu_[regionIdx].appendXPos(Rs); + oilMuTable_[regionIdx].appendXPos(Rs); for (size_t pIdx = 0; pIdx < nP; ++pIdx) { Scalar po = poMin + (poMax - poMin)*pIdx/nP; Scalar muo = muoSpline.eval(po, /*extrapolate=*/true); - oilMu_[regionIdx].appendSamplePoint(RsIdx, po, muo); + oilMuTable_[regionIdx].appendSamplePoint(RsIdx, po, muo); } } } @@ -206,9 +208,9 @@ public: const auto saturatedTable = pvtoTable.getOuterTable(); assert(saturatedTable->numRows() > 1); - auto& oilMu = oilMu_[regionIdx]; - auto& invOilB = inverseOilB_[regionIdx]; - auto& gasDissolutionFactor = gasDissolutionFactor_[regionIdx]; + auto& oilMu = oilMuTable_[regionIdx]; + auto& invOilB = inverseOilBTable_[regionIdx]; + auto& gasDissolutionFactor = gasDissolutionFactorTable_[regionIdx]; gasDissolutionFactor.setXYArrays(saturatedTable->numRows(), saturatedTable->getPressureColumn(), @@ -301,15 +303,15 @@ public: void initEnd() { // calculate the final 2D functions which are used for interpolation. - int numRegions = oilMu_.size(); + int numRegions = oilMuTable_.size(); for (int regionIdx = 0; regionIdx < numRegions; ++ regionIdx) { // calculate the table which stores the inverse of the product of the oil // formation volume factor and the oil viscosity - const auto& oilMu = oilMu_[regionIdx]; - const auto& invOilB = inverseOilB_[regionIdx]; + const auto& oilMu = oilMuTable_[regionIdx]; + const auto& invOilB = inverseOilBTable_[regionIdx]; assert(oilMu.numX() == invOilB.numX()); - auto& invOilBMu = inverseOilBMu_[regionIdx]; + auto& invOilBMu = inverseOilBMuTable_[regionIdx]; for (int rsIdx = 0; rsIdx < oilMu.numX(); ++rsIdx) { invOilBMu.appendXPos(oilMu.xAt(rsIdx)); @@ -328,22 +330,24 @@ public: } } +private: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ - Scalar viscosity(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval viscosity_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { - Scalar Rs = + const LhsEval& Rs = XoG/(1 - XoG) - * BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx) - / BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); + * (BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx) + / BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx)); // ATTENTION: Rs is the first axis! - Scalar invBo = inverseOilB_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true); - Scalar invMuoBo = inverseOilBMu_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true); + const LhsEval& invBo = inverseOilBTable_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true); + const LhsEval& invMuoBo = inverseOilBMuTable_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true); return invBo/invMuoBo; } @@ -351,21 +355,26 @@ public: /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ - Scalar density(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval density_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { Scalar rhooRef = BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx); Scalar rhogRef = BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); + Valgrind::CheckDefined(rhooRef); + Valgrind::CheckDefined(rhogRef); - Scalar Bo = formationVolumeFactor(regionIdx, temperature, pressure, XoG); - Scalar rhoo = rhooRef/Bo; + const LhsEval& Bo = formationVolumeFactor_(regionIdx, temperature, pressure, XoG); + Valgrind::CheckDefined(Bo); + + LhsEval rhoo = rhooRef/Bo; // the oil formation volume factor just represents the partial density of the oil // component in the oil phase. to get the total density of the phase, we have to // add the partial density of the gas component. - Scalar Rs = XoG/(1 - XoG) * rhooRef/rhogRef; + const LhsEval Rs = XoG/(1 - XoG) * rhooRef/rhogRef; rhoo += rhogRef*Rs/Bo; return rhoo; @@ -374,33 +383,37 @@ public: /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ - Scalar formationVolumeFactor(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XoG) const OPM_FINAL + template + LhsEval formationVolumeFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XoG) const { - Scalar Rs = + const LhsEval& Rs = XoG/(1-XoG) - *BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx) - /BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); + *(BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx) + /BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx)); + Valgrind::CheckDefined(Rs); + Valgrind::CheckDefined(XoG); // ATTENTION: Rs is represented by the _first_ axis! - return 1.0 / inverseOilB_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true); + return 1.0 / inverseOilBTable_[regionIdx].eval(Rs, pressure, /*extrapolate=*/true); } /*! * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ - Scalar fugacityCoefficient(int regionIdx, - Scalar temperature, - Scalar pressure, - int compIdx) const OPM_FINAL + template + LhsEval fugacityCoefficient_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int compIdx) const { // set the oil component fugacity coefficient in oil phase // arbitrarily. we use some pseudo-realistic value for the vapor // pressure to ease physical interpretation of the results - Scalar phi_oO = 20e3/pressure; + const LhsEval& phi_oO = 20e3/pressure; if (compIdx == BlackOilFluidSystem::oilCompIdx) return phi_oO; @@ -416,16 +429,15 @@ public: // // first, retrieve the mole fraction of gas a saturated oil // would exhibit at the given pressure - Scalar x_oGSat = saturatedOilGasMoleFraction(regionIdx, temperature, pressure); + const LhsEval& x_oGSat = saturatedOilGasMoleFraction_(regionIdx, temperature, pressure); // then, scale the gas component's gas phase fugacity // coefficient, so that the oil phase ends up at the right // composition if we were doing a flash experiment - Scalar phi_gG = BlackOilFluidSystem::fugCoefficientInGas(gasCompIdx, - temperature, - pressure, - regionIdx); - + const LhsEval& phi_gG = BlackOilFluidSystem::fugCoefficientInGas(gasCompIdx, + temperature, + pressure, + regionIdx); return phi_gG / x_oGSat; } @@ -433,10 +445,11 @@ public: /*! * \brief Returns the gas dissolution factor \f$R_s\f$ [m^3/m^3] of the oil phase. */ - Scalar gasDissolutionFactor(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return gasDissolutionFactor_[regionIdx].eval(pressure, /*extrapolate=*/true); } + template + LhsEval gasDissolutionFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { return gasDissolutionFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true); } /*! * \brief Returns the saturation pressure of the oil phase [Pa] @@ -444,34 +457,38 @@ public: * * \param XoG The mass fraction of the gas component in the oil phase [-] */ - Scalar oilSaturationPressure(int regionIdx, - Scalar temperature, - Scalar XoG) const OPM_FINAL + template + LhsEval oilSaturationPressure_(int regionIdx, + const LhsEval& temperature, + const LhsEval& XoG) const { + typedef Opm::MathToolbox Toolbox; + // use the saturation pressure spline to get a pretty good initial value - Scalar pSat = saturationPressureSpline_[regionIdx].eval(XoG, /*extrapolate=*/true); + LhsEval pSat = saturationPressureSpline_[regionIdx].eval(XoG, /*extrapolate=*/true); + LhsEval eps = pSat*1e-11; // Newton method to do the remaining work. If the initial // value is good, this should only take two to three // iterations... for (int i = 0; i < 20; ++i) { - Scalar f = saturatedOilGasMassFraction(regionIdx, temperature, pSat) - XoG; - Scalar eps = pSat*1e-11; - Scalar fPrime = ((saturatedOilGasMassFraction(regionIdx, temperature, pSat + eps) - XoG) - f)/eps; + const LhsEval& f = saturatedOilGasMassFraction_(regionIdx, temperature, pSat) - XoG; + const LhsEval& fPrime = ((saturatedOilGasMassFraction_(regionIdx, temperature, pSat + eps) - XoG) - f)/eps; - Scalar delta = f/fPrime; + const LhsEval& delta = f/fPrime; pSat -= delta; - if (std::abs(delta) < pSat * 1e-10) + if (std::abs(Toolbox::value(delta)) < Toolbox::value(pSat) * 1e-10) return pSat; } OPM_THROW(NumericalIssue, "Could find the oil saturation pressure for X_o^G = " << XoG); } - Scalar saturatedOilGasMassFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL + template + LhsEval saturatedOilGasMassFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const { Scalar rho_gRef = BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); Scalar rho_oRef = BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx); @@ -479,35 +496,35 @@ public: // calculate the mass of the gas component [kg/m^3] in the oil phase. This is // equivalent to the gas dissolution factor [m^3/m^3] at current pressure times // the gas density [kg/m^3] at standard pressure - Scalar rho_oG = gasDissolutionFactor(regionIdx, temperature, pressure) * rho_gRef; + const LhsEval& rho_oG = gasDissolutionFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true) * rho_gRef; // we now have the total density of saturated oil and the partial density of the // gas component within it. The gas mass fraction is the ratio of these two. return rho_oG/(rho_oRef + rho_oG); } - Scalar saturatedOilGasMoleFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL + template + LhsEval saturatedOilGasMoleFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const { // calculate the mass fractions of gas and oil - Scalar XoG = saturatedOilGasMassFraction(regionIdx, temperature, pressure); + const LhsEval& XoG = saturatedOilGasMassFraction_(regionIdx, temperature, pressure); // which can be converted to mole fractions, given the // components' molar masses Scalar MG = BlackOilFluidSystem::molarMass(gasCompIdx, regionIdx); Scalar MO = BlackOilFluidSystem::molarMass(oilCompIdx, regionIdx); - Scalar avgMolarMass = MO/(1 + XoG*(MO/MG - 1)); - Scalar xoG = XoG*avgMolarMass/MG; + LhsEval avgMolarMass = MO/(1 + XoG*(MO/MG - 1)); + LhsEval xoG = XoG*avgMolarMass/MG; return xoG; } -private: void updateSaturationPressureSpline_(int regionIdx) { - auto& gasDissolutionFactor = gasDissolutionFactor_[regionIdx]; + auto& gasDissolutionFactor = gasDissolutionFactorTable_[regionIdx]; // create the spline representing saturation pressure // depending of the mass fraction in gas @@ -518,7 +535,9 @@ private: Scalar XoG = 0; for (int i=0; i <= n; ++ i) { Scalar pSat = gasDissolutionFactor.xMin() + i*delta; - XoG = saturatedOilGasMassFraction(regionIdx, /*temperature=*/1e100, pSat); + XoG = saturatedOilGasMassFraction_(regionIdx, + /*temperature=*/Scalar(1e100), + pSat); std::pair val(XoG, pSat); pSatSamplePoints.push_back(val); @@ -527,10 +546,10 @@ private: /*type=*/Spline::Monotonic); } - std::vector inverseOilB_; - std::vector oilMu_; - std::vector inverseOilBMu_; - std::vector gasDissolutionFactor_; + std::vector inverseOilBTable_; + std::vector oilMuTable_; + std::vector inverseOilBMuTable_; + std::vector gasDissolutionFactorTable_; std::vector saturationPressureSpline_; }; diff --git a/opm/material/fluidsystems/blackoilpvt/OilPvtInterface.hpp b/opm/material/fluidsystems/blackoilpvt/OilPvtInterface.hpp index c8dfec52f..fbebc866e 100644 --- a/opm/material/fluidsystems/blackoilpvt/OilPvtInterface.hpp +++ b/opm/material/fluidsystems/blackoilpvt/OilPvtInterface.hpp @@ -23,8 +23,9 @@ #ifndef OPM_OIL_PVT_INTERFACE_HPP #define OPM_OIL_PVT_INTERFACE_HPP -namespace Opm { +#include +namespace Opm { /*! * \brief This class represents the Pressure-Volume-Temperature relations of the oil * phase in the black-oil model. @@ -37,9 +38,117 @@ namespace Opm { * Note that, since the application for this class is the black-oil fluid system, the API * exposed by this class is pretty specific to the black-oil model. */ -template +template class OilPvtInterface { +public: + /*! + * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. + */ + virtual Evaluation viscosity(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XoG) const = 0; + virtual Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const = 0; + + /*! + * \brief Returns the formation volume factor [-] of the fluid phase. + */ + virtual Evaluation formationVolumeFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XoG) const = 0; + virtual Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const = 0; + + /*! + * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. + */ + virtual Evaluation density(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XoG) const = 0; + virtual Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const = 0; + + /*! + * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given + * a set of parameters. + */ + virtual Evaluation fugacityCoefficient(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + int compIdx) const = 0; + virtual Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const = 0; + + /*! + * \brief Returns the gas dissolution factor \f$R_s\f$ [m^3/m^3] of saturated oil. + */ + virtual Evaluation gasDissolutionFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; + virtual Scalar gasDissolutionFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; + + /*! + * \brief Returns the saturation pressure [Pa] of oil given the mass fraction of the + * gas component in the oil phase. + * + * Calling this method only makes sense for live oil. All other implementations of + * the black-oil PVT interface will just throw an exception... + */ + virtual Evaluation oilSaturationPressure(int regionIdx, + const Evaluation& temperature, + const Evaluation& XoG) const = 0; + virtual Scalar oilSaturationPressure(int regionIdx, + Scalar temperature, + Scalar XoG) const = 0; + + /*! + * \brief Returns the gas mass fraction of gas-saturated oil at a given temperatire + * and pressure [-]. + * + * Calling this method only makes sense for oil. For all other phases an exception + * will be thrown... + */ + virtual Evaluation saturatedOilGasMassFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; + virtual Scalar saturatedOilGasMassFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; + + /*! + * \brief Returns the gas mole fraction of gas-saturated oil at a given temperatire + * and pressure [-]. + * + * Calling this method only makes sense for oil. For all other phases an exception + * will be thrown... + */ + virtual Evaluation saturatedOilGasMoleFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; + virtual Scalar saturatedOilGasMoleFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; +}; + +// if Evaluation == Scalar, the PVT interface provides only the scalar variants of the +// methods. +template +class OilPvtInterface +{ public: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. @@ -115,6 +224,144 @@ public: Scalar pressure) const = 0; }; +// To prevent the need for most code duplication, derived classes can use this class to +// call templated variants of all the methods +template +class OilPvtInterfaceTemplateWrapper + : public OilPvtInterface +{ + Evaluation viscosity(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XoG) const OPM_FINAL + { return asImp_().template viscosity_(regionIdx, temperature, pressure, XoG); } + Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const OPM_FINAL + { return asImp_().template viscosity_(regionIdx, temperature, pressure, XoG); } + + Evaluation formationVolumeFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XoG) const OPM_FINAL + { return asImp_().template formationVolumeFactor_(regionIdx, temperature, pressure, XoG); } + Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const OPM_FINAL + { return asImp_().template formationVolumeFactor_(regionIdx, temperature, pressure, XoG); } + + Evaluation density(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + const Evaluation& XoG) const OPM_FINAL + { return asImp_().template density_(regionIdx, temperature, pressure, XoG); } + Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const OPM_FINAL + { return asImp_().template density_(regionIdx, temperature, pressure, XoG); } + + Evaluation fugacityCoefficient(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + int compIdx) const OPM_FINAL + { return asImp_().template fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); } + Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const OPM_FINAL + { return asImp_().template fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); } + + Evaluation gasDissolutionFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().template gasDissolutionFactor_(regionIdx, temperature, pressure); } + Scalar gasDissolutionFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().template gasDissolutionFactor_(regionIdx, temperature, pressure); } + + Evaluation oilSaturationPressure(int regionIdx, + const Evaluation& temperature, + const Evaluation& XoG) const OPM_FINAL + { return asImp_().template oilSaturationPressure_(regionIdx, temperature, XoG); } + Scalar oilSaturationPressure(int regionIdx, + Scalar temperature, + Scalar XoG) const OPM_FINAL + { return asImp_().template oilSaturationPressure_(regionIdx, temperature, XoG); } + + Evaluation saturatedOilGasMassFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().template saturatedOilGasMassFraction_(regionIdx, temperature, pressure); } + Scalar saturatedOilGasMassFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().template saturatedOilGasMassFraction_(regionIdx, temperature, pressure); } + + Evaluation saturatedOilGasMoleFraction(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().template saturatedOilGasMoleFraction_(regionIdx, temperature, pressure); } + Scalar saturatedOilGasMoleFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().template saturatedOilGasMoleFraction_(regionIdx, temperature, pressure); } + +private: + const Implementation& asImp_() const + { return *static_cast(this); } +}; + +template +class OilPvtInterfaceTemplateWrapper + : public OilPvtInterface +{ + Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const OPM_FINAL + { return asImp_().template viscosity_(regionIdx, temperature, pressure, XoG); } + Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const OPM_FINAL + { return asImp_().template formationVolumeFactor_(regionIdx, temperature, pressure, XoG); } + Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure, + Scalar XoG) const OPM_FINAL + { return asImp_().template density_(regionIdx, temperature, pressure, XoG); } + Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const OPM_FINAL + { return asImp_().template fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); } + Scalar gasDissolutionFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().template gasDissolutionFactor_(regionIdx, temperature, pressure); } + Scalar oilSaturationPressure(int regionIdx, + Scalar temperature, + Scalar XoG) const OPM_FINAL + { return asImp_().template oilSaturationPressure_(regionIdx, temperature, XoG); } + Scalar saturatedOilGasMassFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().template saturatedOilGasMassFraction_(regionIdx, temperature, pressure); } + Scalar saturatedOilGasMoleFraction(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().template saturatedOilGasMoleFraction_(regionIdx, temperature, pressure); } + +private: + const Implementation& asImp_() const + { return *static_cast(this); } +}; + + } // namespace Opm #endif diff --git a/opm/material/fluidsystems/blackoilpvt/WaterPvtInterface.hpp b/opm/material/fluidsystems/blackoilpvt/WaterPvtInterface.hpp index 58ff95437..da7c1c30c 100644 --- a/opm/material/fluidsystems/blackoilpvt/WaterPvtInterface.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WaterPvtInterface.hpp @@ -23,8 +23,9 @@ #ifndef OPM_WATER_PVT_INTERFACE_HPP #define OPM_WATER_PVT_INTERFACE_HPP -namespace Opm { +#include +namespace Opm { /*! * \brief This class represents the Pressure-Volume-Temperature relations of the water * phase in the black-oil model. @@ -37,13 +38,16 @@ namespace Opm { * Note that, since the application for this class is the black oil fluid system, the API * exposed by this class is pretty specific to the black-oil model. */ -template +template class WaterPvtInterface { public: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ + virtual Evaluation viscosity(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; virtual Scalar viscosity(int regionIdx, Scalar temperature, Scalar pressure) const = 0; @@ -51,6 +55,9 @@ public: /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ + virtual Evaluation formationVolumeFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; virtual Scalar formationVolumeFactor(int regionIdx, Scalar temperature, Scalar pressure) const = 0; @@ -58,6 +65,9 @@ public: /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ + virtual Evaluation density(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const = 0; virtual Scalar density(int regionIdx, Scalar temperature, Scalar pressure) const = 0; @@ -66,12 +76,112 @@ public: * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ + virtual Evaluation fugacityCoefficient(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + int compIdx) const = 0; virtual Scalar fugacityCoefficient(int regionIdx, Scalar temperature, Scalar pressure, int compIdx) const = 0; }; +template +class WaterPvtInterface +{ +public: + virtual Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; + + virtual Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; + virtual Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure) const = 0; + virtual Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const = 0; +}; + +template +class WaterPvtInterfaceTemplateWrapper : public WaterPvtInterface +{ + const Implementation& asImp_() const + { return *static_cast(this); } + +public: + Evaluation viscosity(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().viscosity_(regionIdx, temperature, pressure); } + Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().viscosity_(regionIdx, temperature, pressure); } + + Evaluation formationVolumeFactor(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().formationVolumeFactor_(regionIdx, temperature, pressure); } + Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().formationVolumeFactor_(regionIdx, temperature, pressure); } + + Evaluation density(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure) const OPM_FINAL + { return asImp_().density_(regionIdx, temperature, pressure); } + Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().density_(regionIdx, temperature, pressure); } + + Evaluation fugacityCoefficient(int regionIdx, + const Evaluation& temperature, + const Evaluation& pressure, + int compIdx) const OPM_FINAL + { return asImp_().fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); } + Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const OPM_FINAL + { return asImp_().fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); } +}; + +template +class WaterPvtInterfaceTemplateWrapper + : public WaterPvtInterface +{ + const Implementation& asImp_() const + { return *static_cast(this); } + +public: + Scalar viscosity(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().viscosity_(regionIdx, temperature, pressure); } + + Scalar formationVolumeFactor(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().formationVolumeFactor_(regionIdx, temperature, pressure); } + + Scalar density(int regionIdx, + Scalar temperature, + Scalar pressure) const OPM_FINAL + { return asImp_().density_(regionIdx, temperature, pressure); } + + Scalar fugacityCoefficient(int regionIdx, + Scalar temperature, + Scalar pressure, + int compIdx) const OPM_FINAL + { return asImp_().fugacityCoefficient_(regionIdx, temperature, pressure, compIdx); } +}; + } // namespace Opm #endif diff --git a/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp b/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp index ee4656ec1..c09059b37 100644 --- a/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp +++ b/opm/material/fluidsystems/blackoilpvt/WetGasPvt.hpp @@ -35,15 +35,17 @@ #endif namespace Opm { - /*! * \brief This class represents the Pressure-Volume-Temperature relations of the gas phas * with vaporized oil. */ -template -class WetGasPvt : public GasPvtInterface +template +class WetGasPvt + : public GasPvtInterfaceTemplateWrapper > { - typedef FluidSystems::BlackOil BlackOilFluidSystem; + friend class GasPvtInterfaceTemplateWrapper >; + + typedef FluidSystems::BlackOil BlackOilFluidSystem; typedef Opm::UniformXTabulated2DFunction TabulatedTwoDFunction; typedef Opm::Tabulated1DFunction TabulatedOneDFunction; @@ -64,7 +66,7 @@ public: inverseGasB_.resize(numRegions); inverseGasBMu_.resize(numRegions); gasMu_.resize(numRegions); - oilVaporizationFactor_.resize(numRegions); + oilVaporizationFactorTable_.resize(numRegions); saturationPressureSpline_.resize(numRegions); } @@ -74,7 +76,7 @@ public: * \param samplePoints A container of (x,y) values. */ void setSaturatedGasOilVaporizationFactor(int regionIdx, const SamplingPoints &samplePoints) - { oilVaporizationFactor_[regionIdx].setContainerOfTuples(samplePoints); } + { oilVaporizationFactorTable_[regionIdx].setContainerOfTuples(samplePoints); } /*! * \brief Initialize the function for the gas formation volume factor @@ -89,12 +91,12 @@ public: { auto& invGasB = inverseGasB_[regionIdx]; - auto &Rv = oilVaporizationFactor_[regionIdx]; + auto &Rv = oilVaporizationFactorTable_[regionIdx]; Scalar T = BlackOilFluidSystem::surfaceTemperature; Scalar RvMin = 0.0; - Scalar RvMax = Rv.eval(oilVaporizationFactor_[regionIdx].xMax(), /*extrapolate=*/true); + Scalar RvMax = Rv.eval(oilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true); Scalar poMin = samplePoints.front().first; Scalar poMax = samplePoints.back().first; @@ -165,10 +167,10 @@ public: */ void setSaturatedGasViscosity(int regionIdx, const SamplingPoints &samplePoints ) { - auto& oilVaporizationFactor = oilVaporizationFactor_[regionIdx]; + auto& oilVaporizationFactor = oilVaporizationFactorTable_[regionIdx]; Scalar RvMin = 0.0; - Scalar RvMax = oilVaporizationFactor.eval(oilVaporizationFactor_[regionIdx].xMax(), /*extrapolate=*/true); + Scalar RvMax = oilVaporizationFactor.eval(oilVaporizationFactorTable_[regionIdx].xMax(), /*extrapolate=*/true); Scalar poMin = samplePoints.front().first; Scalar poMax = samplePoints.back().first; @@ -206,11 +208,11 @@ public: auto& gasMu = gasMu_[regionIdx]; auto& invGasB = inverseGasB_[regionIdx]; - auto& oilVaporizationFactor = oilVaporizationFactor_[regionIdx]; + auto& oilVaporizationFactor = oilVaporizationFactorTable_[regionIdx]; oilVaporizationFactor.setXYArrays(saturatedTable->numRows(), - saturatedTable->getPressureColumn(), - saturatedTable->getOilSolubilityColumn()); + saturatedTable->getPressureColumn(), + saturatedTable->getOilSolubilityColumn()); // extract the table for the gas dissolution and the oil formation volume factors for (int outerIdx = 0; outerIdx < static_cast(saturatedTable->numRows()); ++ outerIdx) { @@ -326,21 +328,23 @@ public: } } +private: /*! * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters. */ - Scalar viscosity(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XgO) const OPM_FINAL + template + LhsEval viscosity_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO) const { - Scalar Rv = + const LhsEval& Rv = XgO/(1 - XgO) - * BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx) - / BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); + * (BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx) + / BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx)); - Scalar invBg = inverseGasB_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true); - Scalar invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true); + const LhsEval& invBg = inverseGasB_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true); + const LhsEval& invMugBg = inverseGasBMu_[regionIdx].eval(pressure, Rv, /*extrapolate=*/true); return invBg/invMugBg; } @@ -348,22 +352,23 @@ public: /*! * \brief Returns the density [kg/m^3] of the fluid phase given a set of parameters. */ - Scalar density(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XgO) const OPM_FINAL + template + LhsEval density_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO) const { Scalar rhooRef = BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx); Scalar rhogRef = BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); - Scalar Bg = formationVolumeFactor(regionIdx, temperature, pressure, XgO); - Scalar rhoo = rhooRef/Bg; + const LhsEval& Bg = formationVolumeFactor_(regionIdx, temperature, pressure, XgO); + LhsEval rhoo = rhooRef/Bg; // the oil formation volume factor just represents the partial density of the gas // component in the gas phase. to get the total density of the phase, we have to // add the partial density of the oil component. - Scalar Rv = XgO/(1 - XgO) * rhooRef/rhogRef; - rhoo += rhogRef*Rv/Bg; + const LhsEval& Rv = XgO/(1 - XgO) * (rhooRef/rhogRef); + rhoo += (rhogRef*Rv)/Bg; return rhoo; } @@ -371,12 +376,13 @@ public: /*! * \brief Returns the formation volume factor [-] of the fluid phase. */ - Scalar formationVolumeFactor(int regionIdx, - Scalar temperature, - Scalar pressure, - Scalar XgO) const OPM_FINAL + template + LhsEval formationVolumeFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + const LhsEval& XgO) const { - Scalar Rv = + const LhsEval& Rv = XgO/(1-XgO) *BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx) /BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); @@ -388,15 +394,18 @@ public: * \brief Returns the fugacity coefficient [Pa] of a component in the fluid phase given * a set of parameters. */ - Scalar fugacityCoefficient(int regionIdx, - Scalar temperature, - Scalar pressure, - int compIdx) const OPM_FINAL + template + LhsEval fugacityCoefficient_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure, + int compIdx) const { + typedef Opm::MathToolbox Toolbox; + // set the oil component fugacity coefficient in oil phase // arbitrarily. we use some pseudo-realistic value for the vapor // pressure to ease physical interpretation of the results - Scalar phi_gG = 1.0; + LhsEval phi_gG = Toolbox::createConstant(1.0); if (compIdx == BlackOilFluidSystem::gasCompIdx) return phi_gG; @@ -411,16 +420,15 @@ public: // // first, retrieve the mole fraction of gas a saturated oil // would exhibit at the given pressure - Scalar x_gOSat = saturatedGasOilMoleFraction(regionIdx, temperature, pressure); + const LhsEval& x_gOSat = saturatedGasOilMoleFraction_(regionIdx, temperature, pressure); // then, scale the oil component's gas phase fugacity // coefficient, so that the oil phase ends up at the right // composition if we were doing a flash experiment - Scalar phi_oO = BlackOilFluidSystem::fugCoefficientInOil(oilCompIdx, - temperature, - pressure, - regionIdx); - + const LhsEval& phi_oO = BlackOilFluidSystem::fugCoefficientInOil(oilCompIdx, + temperature, + pressure, + regionIdx); return phi_oO / x_gOSat; } @@ -428,10 +436,11 @@ public: /*! * \brief Returns the gas dissolution factor \f$R_s\f$ [m^3/m^3] of the oil phase. */ - Scalar oilVaporizationFactor(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL - { return oilVaporizationFactor_[regionIdx].eval(pressure, /*extrapolate=*/true); } + template + LhsEval oilVaporizationFactor_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const + { return oilVaporizationFactorTable_[regionIdx].eval(pressure, /*extrapolate=*/true); } /*! * \brief Returns the saturation pressure of the gas phase [Pa] @@ -439,34 +448,38 @@ public: * * \param XgO The mass fraction of the oil component in the gas phase [-] */ - Scalar gasSaturationPressure(int regionIdx, - Scalar temperature, - Scalar XgO) const OPM_FINAL + template + LhsEval gasSaturationPressure_(int regionIdx, + const LhsEval& temperature, + const LhsEval& XgO) const { + typedef Opm::MathToolbox Toolbox; + // use the saturation pressure spline to get a pretty good initial value - Scalar pSat = saturationPressureSpline_[regionIdx].eval(XgO, /*extrapolate=*/true); + LhsEval pSat = saturationPressureSpline_[regionIdx].eval(XgO, /*extrapolate=*/true); + const LhsEval& eps = pSat*1e-11; // Newton method to do the remaining work. If the initial // value is good, this should only take two to three // iterations... for (int i = 0; i < 20; ++i) { - Scalar f = saturatedGasOilMassFraction(regionIdx, temperature, pSat) - XgO; - Scalar eps = pSat*1e-11; - Scalar fPrime = ((saturatedGasOilMassFraction(regionIdx, temperature, pSat + eps) - XgO) - f)/eps; + const LhsEval& f = saturatedGasOilMassFraction_(regionIdx, temperature, pSat) - XgO; + const LhsEval& fPrime = ((saturatedGasOilMassFraction_(regionIdx, temperature, pSat + eps) - XgO) - f)/eps; - Scalar delta = f/fPrime; + const LhsEval& delta = f/fPrime; pSat -= delta; - if (std::abs(delta) < pSat * 1e-10) + if (std::abs(Toolbox::value(delta)) < std::abs(Toolbox::value(pSat)) * 1e-10) return pSat; } OPM_THROW(NumericalIssue, "Could find the gas saturation pressure for X_g^O = " << XgO); } - Scalar saturatedGasOilMassFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL + template + LhsEval saturatedGasOilMassFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const { Scalar rho_gRef = BlackOilFluidSystem::referenceDensity(gasPhaseIdx, regionIdx); Scalar rho_oRef = BlackOilFluidSystem::referenceDensity(oilPhaseIdx, regionIdx); @@ -474,27 +487,28 @@ public: // calculate the mass of the oil component [kg/m^3] in the oil phase. This is // equivalent to the gas dissolution factor [m^3/m^3] at current pressure times // the gas density [kg/m^3] at standard pressure - Scalar rho_oG = oilVaporizationFactor(regionIdx, temperature, pressure) * rho_gRef; + const LhsEval& rho_oG = oilVaporizationFactor_(regionIdx, temperature, pressure) * rho_gRef; // we now have the total density of saturated oil and the partial density of the // oil component within it. The gas mass fraction is the ratio of these two. return rho_oG/(rho_oRef + rho_oG); } - Scalar saturatedGasOilMoleFraction(int regionIdx, - Scalar temperature, - Scalar pressure) const OPM_FINAL + template + LhsEval saturatedGasOilMoleFraction_(int regionIdx, + const LhsEval& temperature, + const LhsEval& pressure) const { // calculate the mass fractions of gas and oil - Scalar XgO = saturatedGasOilMassFraction(regionIdx, temperature, pressure); + const LhsEval& XgO = saturatedGasOilMassFraction_(regionIdx, temperature, pressure); // which can be converted to mole fractions, given the // components' molar masses Scalar MG = BlackOilFluidSystem::molarMass(gasCompIdx, regionIdx); Scalar MO = BlackOilFluidSystem::molarMass(oilCompIdx, regionIdx); - Scalar avgMolarMass = MO/(1 + (1 - XgO)*(MO/MG - 1)); - Scalar xgO = XgO*avgMolarMass/MG; + const LhsEval& avgMolarMass = MO/(1 + (1 - XgO)*(MO/MG - 1)); + const LhsEval& xgO = XgO*avgMolarMass/MG; return xgO; } @@ -502,7 +516,7 @@ public: private: void updateSaturationPressureSpline_(int regionIdx) { - auto& oilVaporizationFactor = oilVaporizationFactor_[regionIdx]; + auto& oilVaporizationFactor = oilVaporizationFactorTable_[regionIdx]; // create the spline representing saturation pressure // depending of the mass fraction in gas @@ -513,7 +527,7 @@ private: Scalar XgO = 0; for (int i=0; i <= n; ++ i) { Scalar pSat = oilVaporizationFactor.xMin() + i*delta; - XgO = saturatedGasOilMassFraction(regionIdx, /*temperature=*/1e100, pSat); + XgO = saturatedGasOilMassFraction_(regionIdx, /*temperature=*/Scalar(1e100), pSat); std::pair val(XgO, pSat); pSatSamplePoints.push_back(val); @@ -525,7 +539,7 @@ private: std::vector inverseGasB_; std::vector gasMu_; std::vector inverseGasBMu_; - std::vector oilVaporizationFactor_; + std::vector oilVaporizationFactorTable_; std::vector saturationPressureSpline_; }; diff --git a/tests/checkFluidSystem.hpp b/tests/checkFluidSystem.hpp index 7697072ec..1e3305966 100644 --- a/tests/checkFluidSystem.hpp +++ b/tests/checkFluidSystem.hpp @@ -253,7 +253,7 @@ void checkFluidState(const BaseFluidState &fs) /*! * \brief Checks whether a fluid system adheres to the specification. */ -template +template void checkFluidSystem() { std::cout << "Testing fluid system '" << Opm::className() << "'\n"; @@ -263,7 +263,8 @@ void checkFluidSystem() enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; - HairSplittingFluidState fs; + typedef HairSplittingFluidState FluidState; + FluidState fs; fs.allowTemperature(true); fs.allowPressure(true); fs.allowComposition(true); @@ -290,7 +291,8 @@ void checkFluidSystem() // some value to make sure the return values of the fluid system // are convertible to scalars - Scalar OPM_UNUSED val; + LhsEval OPM_UNUSED val; + Scalar OPM_UNUSED scalarVal; // actually check the fluid system API try { FluidSystem::init(); } catch (...) {}; @@ -299,22 +301,35 @@ void checkFluidSystem() fs.allowPressure(FluidSystem::isCompressible(phaseIdx)); fs.allowComposition(true); fs.allowDensity(false); - try { val = FluidSystem::density(fs, paramCache, phaseIdx); } catch (...) {}; + try { auto OPM_UNUSED tmpVal = FluidSystem::density(fs, paramCache, phaseIdx); static_assert(std::is_same::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {}; + try { val = FluidSystem::template density(fs, paramCache, phaseIdx); } catch (...) {}; + try { scalarVal = FluidSystem::template density(fs, paramCache, phaseIdx); } catch (...) {}; fs.allowPressure(true); fs.allowDensity(true); - try { val = FluidSystem::viscosity(fs, paramCache, phaseIdx); } catch (...) {}; - try { val = FluidSystem::enthalpy(fs, paramCache, phaseIdx); } catch (...) {}; - try { val = FluidSystem::heatCapacity(fs, paramCache, phaseIdx); } catch (...) {}; - try { val = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); } catch (...) {}; + try { auto OPM_UNUSED tmpVal = FluidSystem::viscosity(fs, paramCache, phaseIdx); static_assert(std::is_same::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {}; + try { auto OPM_UNUSED tmpVal = FluidSystem::enthalpy(fs, paramCache, phaseIdx); static_assert(std::is_same::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {}; + try { auto OPM_UNUSED tmpVal = FluidSystem::heatCapacity(fs, paramCache, phaseIdx); static_assert(std::is_same::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {}; + try { auto OPM_UNUSED tmpVal = FluidSystem::thermalConductivity(fs, paramCache, phaseIdx); static_assert(std::is_same::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {}; + try { val = FluidSystem::template viscosity(fs, paramCache, phaseIdx); } catch (...) {}; + try { val = FluidSystem::template enthalpy(fs, paramCache, phaseIdx); } catch (...) {}; + try { val = FluidSystem::template heatCapacity(fs, paramCache, phaseIdx); } catch (...) {}; + try { val = FluidSystem::template thermalConductivity(fs, paramCache, phaseIdx); } catch (...) {}; + try { scalarVal = FluidSystem::template viscosity(fs, paramCache, phaseIdx); } catch (...) {}; + try { scalarVal = FluidSystem::template enthalpy(fs, paramCache, phaseIdx); } catch (...) {}; + try { scalarVal = FluidSystem::template heatCapacity(fs, paramCache, phaseIdx); } catch (...) {}; + try { scalarVal = FluidSystem::template thermalConductivity(fs, paramCache, phaseIdx); } catch (...) {}; for (int compIdx = 0; compIdx < numComponents; ++ compIdx) { fs.allowComposition(!FluidSystem::isIdealMixture(phaseIdx)); - try { val = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; + try { auto OPM_UNUSED tmpVal = FluidSystem::fugacityCoefficient(fs, paramCache, phaseIdx, compIdx); static_assert(std::is_same::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {}; + try { val = FluidSystem::template fugacityCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; + try { scalarVal = FluidSystem::template fugacityCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; fs.allowComposition(true); - try { val = FluidSystem::diffusionCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; + try { auto OPM_UNUSED tmpVal = FluidSystem::diffusionCoefficient(fs, paramCache, phaseIdx, compIdx); static_assert(std::is_same::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {}; + try { val = FluidSystem::template diffusionCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; + try { scalarVal = FluidSystem::template fugacityCoefficient(fs, paramCache, phaseIdx, compIdx); } catch (...) {}; } - } // test for phaseName(), isLiquid() and isIdealGas() @@ -324,7 +339,7 @@ void checkFluidSystem() bVal = FluidSystem::isIdealGas(phaseIdx); } - // test for componentName() + // test for molarMass() and componentName() for (int compIdx = 0; compIdx < numComponents; ++ compIdx) { val = FluidSystem::molarMass(compIdx); std::string OPM_UNUSED name = FluidSystem::componentName(compIdx); diff --git a/tests/test_fluidsystems.cpp b/tests/test_fluidsystems.cpp index ae381de5b..c58e826ca 100644 --- a/tests/test_fluidsystems.cpp +++ b/tests/test_fluidsystems.cpp @@ -24,6 +24,9 @@ */ #include "config.h" +#include +#include + #include "checkFluidSystem.hpp" // include all fluid systems in opm-material @@ -45,8 +48,7 @@ #include #include -// include the tables for CO2 which are delivered with opm-material by -// default +// include the tables for CO2 which are delivered with opm-material by default #include namespace Opm { @@ -77,116 +79,137 @@ public: }; }; -int main(int argc, char **argv) +// check the API of all fluid states +template +void testAllFluidStates() +{ + typedef Opm::FluidSystems::H2ON2 FluidSystem; + + // CompositionalFluidState + { Opm::CompositionalFluidState fs; + checkFluidState(fs); } + + // NonEquilibriumFluidState + { Opm::NonEquilibriumFluidState fs; + checkFluidState(fs); } + + // ImmiscibleFluidState + { Opm::ImmiscibleFluidState fs; + checkFluidState(fs); } + + typedef Opm::CompositionalFluidState BaseFluidState; + BaseFluidState baseFs; + + // TemperatureOverlayFluidState + { Opm::TemperatureOverlayFluidState fs(baseFs); + checkFluidState(fs); } + + // PressureOverlayFluidState + { Opm::PressureOverlayFluidState fs(baseFs); + checkFluidState(fs); } + + // SaturationOverlayFluidState + { Opm::SaturationOverlayFluidState fs(baseFs); + checkFluidState(fs); } +} + +template +void testAllFluidSystems() { - typedef double Scalar; typedef Opm::H2O H2O; typedef Opm::N2 N2; typedef Opm::LiquidPhase Liquid; typedef Opm::GasPhase Gas; - MyMpiHelper mpiHelper(argc, argv); - - // check all fluid states - { - typedef Opm::FluidSystems::H2ON2 FluidSystem; - - // CompositionalFluidState - { Opm::CompositionalFluidState fs; - checkFluidState(fs); } - - // NonEquilibriumFluidState - { Opm::NonEquilibriumFluidState fs; - checkFluidState(fs); } - - // ImmiscibleFluidState - { Opm::ImmiscibleFluidState fs; - checkFluidState(fs); } - - typedef Opm::CompositionalFluidState BaseFluidState; - BaseFluidState baseFs; - - // TemperatureOverlayFluidState - { Opm::TemperatureOverlayFluidState fs(baseFs); - checkFluidState(fs); } - - // PressureOverlayFluidState - { Opm::PressureOverlayFluidState fs(baseFs); - checkFluidState(fs); } - - // SaturationOverlayFluidState - { Opm::SaturationOverlayFluidState fs(baseFs); - checkFluidState(fs); } - } - // black-oil - { typedef Opm::FluidSystems::BlackOil FluidSystem; - if (false) checkFluidSystem(); } + { typedef Opm::FluidSystems::BlackOil FluidSystem; + if (false) checkFluidSystem(); } // Brine -- CO2 { typedef Opm::FluidSystems::BrineCO2 FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } // H2O -- N2 { typedef Opm::FluidSystems::H2ON2 FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::FluidSystems::H2ON2 FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } // H2O -- N2 -- liquid phase { typedef Opm::FluidSystems::H2ON2LiquidPhase FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::FluidSystems::H2ON2LiquidPhase FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } // H2O -- Air { typedef Opm::SimpleH2O H2O; const bool enableComplexRelations=false; typedef Opm::FluidSystems::H2OAir FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::SimpleH2O H2O; const bool enableComplexRelations=true; typedef Opm::FluidSystems::H2OAir FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::H2O H2O; const bool enableComplexRelations=false; typedef Opm::FluidSystems::H2OAir FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::H2O H2O; const bool enableComplexRelations=true; typedef Opm::FluidSystems::H2OAir FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } // H2O -- Air -- Mesitylene { typedef Opm::FluidSystems::H2OAirMesitylene FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } // H2O -- Air -- Xylene { typedef Opm::FluidSystems::H2OAirXylene FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } // 2p-immiscible { typedef Opm::FluidSystems::TwoPhaseImmiscible FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::FluidSystems::TwoPhaseImmiscible FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::FluidSystems::TwoPhaseImmiscible FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } // 1p { typedef Opm::FluidSystems::SinglePhase FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } { typedef Opm::FluidSystems::SinglePhase FluidSystem; - checkFluidSystem(); } + checkFluidSystem(); } +} + +class TestAdTag; + +int main(int argc, char **argv) +{ + typedef double Scalar; + typedef Opm::LocalAd::Evaluation Evaluation; + + MyMpiHelper mpiHelper(argc, argv); + + // ensure that all fluid states are API-compliant + testAllFluidStates(); + testAllFluidStates(); + + // ensure that all fluid systems are API-compliant: Each fluid system must be usable + // for both, scalars and function evaluations. The fluid systems for function + // evaluations must also be usable for scalars. + testAllFluidSystems(); + testAllFluidSystems(); + testAllFluidSystems(); return 0; }