From e009e640c5c42af58dd4a08a04e68bd4ba2ae693 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 21 May 2015 15:33:14 +0200 Subject: [PATCH] make the fluid states ready for the local-AD framework in fact they wouldn't have needed any modification, but returning constant references instead of copies of the stored values saves quite a few calls to copy constructors. besides this, fluid states are now required to export the 'Scalar' type, which allows to use it outside of the FluidState without resorting to the c++ 'decltype' construct. --- .../fluidmatrixinteractions/EffToAbsLaw.hpp | 12 +++---- .../FluidStateCompositionModules.hpp | 24 ++++++++----- .../fluidstates/FluidStateDensityModules.hpp | 15 ++++---- .../fluidstates/FluidStateEnthalpyModules.hpp | 35 ++++++++++++------- .../fluidstates/FluidStateFugacityModules.hpp | 10 +++--- .../fluidstates/FluidStatePressureModules.hpp | 12 ++++--- .../FluidStateSaturationModules.hpp | 15 ++++---- .../FluidStateTemperatureModules.hpp | 9 +++-- .../FluidStateViscosityModules.hpp | 10 +++--- .../fluidstates/PressureOverlayFluidState.hpp | 4 ++- .../SaturationOverlayFluidState.hpp | 4 ++- .../TemperatureOverlayFluidState.hpp | 4 ++- tests/checkFluidSystem.hpp | 11 ++++-- tests/test_fluidsystems.cpp | 6 ++-- 14 files changed, 109 insertions(+), 62 deletions(-) diff --git a/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp b/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp index e2b50fe5c..ac2998209 100644 --- a/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp +++ b/opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp @@ -112,7 +112,7 @@ public: template static void capillaryPressures(Container &values, const Params ¶ms, const FluidState &fs) { - typedef Opm::SaturationOverlayFluidState OverlayFluidState; + typedef Opm::SaturationOverlayFluidState OverlayFluidState; OverlayFluidState overlayFs(fs); for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { @@ -138,7 +138,7 @@ public: template static void relativePermeabilities(Container &values, const Params ¶ms, const FluidState &fs) { - typedef Opm::SaturationOverlayFluidState OverlayFluidState; + typedef Opm::SaturationOverlayFluidState OverlayFluidState; OverlayFluidState overlayFs(fs); for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) { @@ -165,7 +165,7 @@ public: template static Scalar pcnw(const Params ¶ms, const FluidState &fs) { - typedef Opm::SaturationOverlayFluidState OverlayFluidState; + typedef Opm::SaturationOverlayFluidState OverlayFluidState; static_assert(FluidState::numPhases == numPhases, "The fluid state and the material law must exhibit the same " @@ -294,7 +294,7 @@ public: template static Scalar krw(const Params ¶ms, const FluidState &fs) { - typedef Opm::SaturationOverlayFluidState OverlayFluidState; + typedef Opm::SaturationOverlayFluidState OverlayFluidState; static_assert(FluidState::numPhases == numPhases, "The fluid state and the material law must exhibit the same " @@ -322,7 +322,7 @@ public: template static Scalar krn(const Params ¶ms, const FluidState &fs) { - typedef Opm::SaturationOverlayFluidState OverlayFluidState; + typedef Opm::SaturationOverlayFluidState OverlayFluidState; static_assert(FluidState::numPhases == numPhases, "The fluid state and the material law must exhibit the same " @@ -353,7 +353,7 @@ public: static typename std::enable_if< (Traits::numPhases > 2), ScalarT>::type krg(const Params ¶ms, const FluidState &fs) { - typedef Opm::SaturationOverlayFluidState OverlayFluidState; + typedef Opm::SaturationOverlayFluidState OverlayFluidState; static_assert(FluidState::numPhases == numPhases, "The fluid state and the material law must exhibit the same " diff --git a/opm/material/fluidstates/FluidStateCompositionModules.hpp b/opm/material/fluidstates/FluidStateCompositionModules.hpp index eba0c6787..988f900ef 100644 --- a/opm/material/fluidstates/FluidStateCompositionModules.hpp +++ b/opm/material/fluidstates/FluidStateCompositionModules.hpp @@ -25,6 +25,7 @@ #define OPM_FLUID_STATE_COMPOSITION_MODULES_HPP #include +#include #include #include @@ -57,7 +58,7 @@ public: /*! * \brief The mole fraction of a component in a phase [] */ - Scalar moleFraction(int phaseIdx, int compIdx) const + const Scalar& moleFraction(int phaseIdx, int compIdx) const { return moleFraction_[phaseIdx][compIdx]; } /*! @@ -65,11 +66,13 @@ public: */ Scalar massFraction(int phaseIdx, int compIdx) const { + typedef Opm::MathToolbox Toolbox; + return - std::abs(sumMoleFractions_[phaseIdx]) - * moleFraction_[phaseIdx][compIdx] - * FluidSystem::molarMass(compIdx) - / std::max(1e-40, std::abs(averageMolarMass_[phaseIdx])); + Toolbox::abs(sumMoleFractions_[phaseIdx]) + *moleFraction_[phaseIdx][compIdx] + *FluidSystem::molarMass(compIdx) + / Toolbox::max(1e-40, Toolbox::abs(averageMolarMass_[phaseIdx])); } /*! @@ -80,7 +83,7 @@ public: * component's molar masses weighted by the current mole fraction: * \f[ \bar M_\alpha = \sum_\kappa M^\kappa x_\alpha^\kappa \f] */ - Scalar averageMolarMass(int phaseIdx) const + const Scalar& averageMolarMass(int phaseIdx) const { return averageMolarMass_[phaseIdx]; } /*! @@ -100,7 +103,7 @@ public: * and update the average molar mass [kg/mol] according * to the current composition of the phase */ - void setMoleFraction(int phaseIdx, int compIdx, Scalar value) + void setMoleFraction(int phaseIdx, int compIdx, const Scalar& value) { Valgrind::CheckDefined(value); Valgrind::SetUndefined(sumMoleFractions_[phaseIdx]); @@ -125,11 +128,16 @@ public: template void assign(const FluidState& fs) { + typedef typename FluidState::Scalar FsScalar; + typedef Opm::MathToolbox FsToolbox; + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { averageMolarMass_[phaseIdx] = 0; sumMoleFractions_[phaseIdx] = 0; for (int compIdx = 0; compIdx < numComponents; ++compIdx) { - moleFraction_[phaseIdx][compIdx] = fs.moleFraction(phaseIdx, compIdx); + moleFraction_[phaseIdx][compIdx] = + FsToolbox::template toLhs(fs.moleFraction(phaseIdx, compIdx)); + averageMolarMass_[phaseIdx] += moleFraction_[phaseIdx][compIdx]*FluidSystem::molarMass(compIdx); sumMoleFractions_[phaseIdx] += moleFraction_[phaseIdx][compIdx]; } diff --git a/opm/material/fluidstates/FluidStateDensityModules.hpp b/opm/material/fluidstates/FluidStateDensityModules.hpp index abc8e089d..f116c0e29 100644 --- a/opm/material/fluidstates/FluidStateDensityModules.hpp +++ b/opm/material/fluidstates/FluidStateDensityModules.hpp @@ -24,11 +24,12 @@ #ifndef OPM_FLUID_STATE_DENSITY_MODULES_HPP #define OPM_FLUID_STATE_DENSITY_MODULES_HPP -#include - #include #include +#include +#include + #include namespace Opm { @@ -51,7 +52,7 @@ public: /*! * \brief The density of a fluid phase [kg/m^3] */ - Scalar density(int phaseIdx) const + const Scalar& density(int phaseIdx) const { return density_[phaseIdx]; } /*! @@ -69,7 +70,7 @@ public: /*! * \brief Set the density of a phase [kg/m^3] */ - void setDensity(int phaseIdx, Scalar value) + void setDensity(int phaseIdx, const Scalar& value) { density_[phaseIdx] = value; } /*! @@ -79,8 +80,10 @@ public: template void assign(const FluidState& fs) { + typedef typename FluidState::Scalar FsScalar; + typedef Opm::MathToolbox FsToolbox; for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - density_[phaseIdx] = fs.density(phaseIdx); + density_[phaseIdx] = FsToolbox::template toLhs(fs.density(phaseIdx)); } } @@ -120,7 +123,7 @@ public: /*! * \brief The density of a fluid phase [kg/m^3] */ - Scalar density(int phaseIdx) const + const Scalar& density(int phaseIdx) const { OPM_THROW(std::logic_error, "Density is not provided by this fluid state"); } diff --git a/opm/material/fluidstates/FluidStateEnthalpyModules.hpp b/opm/material/fluidstates/FluidStateEnthalpyModules.hpp index fbadcc817..736dca8d4 100644 --- a/opm/material/fluidstates/FluidStateEnthalpyModules.hpp +++ b/opm/material/fluidstates/FluidStateEnthalpyModules.hpp @@ -24,15 +24,15 @@ #ifndef OPM_FLUID_STATE_ENTHALPY_MODULES_HPP #define OPM_FLUID_STATE_ENTHALPY_MODULES_HPP -#include - #include #include +#include +#include + #include namespace Opm { - /*! * \brief Module for the modular fluid state which stores the * enthalpies explicitly. @@ -51,7 +51,7 @@ public: /*! * \brief The specific enthalpy of a fluid phase [J/kg] */ - Scalar enthalpy(int phaseIdx) const + const Scalar& enthalpy(int phaseIdx) const { return enthalpy_[phaseIdx]; } /*! @@ -63,7 +63,7 @@ public: /*! * \brief Set the specific enthalpy of a phase [J/kg] */ - void setEnthalpy(int phaseIdx, Scalar value) + void setEnthalpy(int phaseIdx, const Scalar& value) { enthalpy_[phaseIdx] = value; } /*! @@ -73,8 +73,10 @@ public: template void assign(const FluidState& fs) { + typedef typename FluidState::Scalar FsScalar; + typedef Opm::MathToolbox FsToolbox; for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - enthalpy_[phaseIdx] = fs.enthalpy(phaseIdx); + enthalpy_[phaseIdx] = FsToolbox::template toLhs(fs.enthalpy(phaseIdx)); } } @@ -100,8 +102,9 @@ protected: /*! * \brief Module for the modular fluid state which does not store the - * enthalpies but throws std::logic_error instead. - */ + * enthalpies but returns 0 instead. + * + * Also, the returned values are marked as undefined in Valgrind... */ template @@ -114,14 +117,22 @@ public: /*! * \brief The specific internal energy of a fluid phase [J/kg] */ - Scalar internalEnergy(int phaseIdx) const - { OPM_THROW(std::logic_error, "Internal energy is not provided by this fluid state"); } + const Scalar& internalEnergy(int phaseIdx) const + { + static Scalar tmp = 0; + Valgrind::SetUndefined(tmp); + return tmp; + } /*! * \brief The specific enthalpy of a fluid phase [J/kg] */ - Scalar enthalpy(int phaseIdx) const - { OPM_THROW(std::logic_error, "Enthalpy is not provided by this fluid state"); } + const Scalar& enthalpy(int phaseIdx) const + { + static Scalar tmp = 0; + Valgrind::SetUndefined(tmp); + return tmp; + } /*! * \brief Retrieve all parameters from an arbitrary fluid diff --git a/opm/material/fluidstates/FluidStateFugacityModules.hpp b/opm/material/fluidstates/FluidStateFugacityModules.hpp index fa82f2bd4..97b871c3d 100644 --- a/opm/material/fluidstates/FluidStateFugacityModules.hpp +++ b/opm/material/fluidstates/FluidStateFugacityModules.hpp @@ -55,7 +55,7 @@ public: /*! * \brief The fugacity coefficient of a component in a phase [] */ - Scalar fugacityCoefficient(int phaseIdx, int compIdx) const + const Scalar& fugacityCoefficient(int phaseIdx, int compIdx) const { return fugacityCoefficient_[phaseIdx][compIdx]; } /*! @@ -67,7 +67,7 @@ public: /*! * \brief Set the fugacity of a component in a phase [] */ - void setFugacityCoefficient(int phaseIdx, int compIdx, Scalar value) + void setFugacityCoefficient(int phaseIdx, int compIdx, const Scalar& value) { fugacityCoefficient_[phaseIdx][compIdx] = value; } /*! @@ -137,7 +137,7 @@ public: /*! * \brief Set the fugacity of a component in a phase [] */ - void setFugacityCoefficient(int phaseIdx, Scalar value) + void setFugacityCoefficient(int phaseIdx, const Scalar& value) { fugacityCoefficient_[phaseIdx] = value; } /*! @@ -188,13 +188,13 @@ public: /*! * \brief The fugacity coefficient of a component in a phase [] */ - Scalar fugacityCoefficient(int phaseIdx, int compIdx) const + const Scalar& fugacityCoefficient(int phaseIdx, int compIdx) const { OPM_THROW(std::logic_error, "Fugacity coefficients are not provided by this fluid state"); } /*! * \brief The fugacity of a component in a phase [Pa] */ - Scalar fugacity(int phaseIdx, int compIdx) const + const Scalar& fugacity(int phaseIdx, int compIdx) const { OPM_THROW(std::logic_error, "Fugacities coefficients are not provided by this fluid state"); } /*! diff --git a/opm/material/fluidstates/FluidStatePressureModules.hpp b/opm/material/fluidstates/FluidStatePressureModules.hpp index 723aa33d9..9fff6db2f 100644 --- a/opm/material/fluidstates/FluidStatePressureModules.hpp +++ b/opm/material/fluidstates/FluidStatePressureModules.hpp @@ -24,8 +24,8 @@ #ifndef OPM_FLUID_STATE_PRESSURE_MODULES_HPP #define OPM_FLUID_STATE_PRESSURE_MODULES_HPP +#include #include - #include #include @@ -51,14 +51,14 @@ public: /*! * \brief The pressure of a fluid phase [Pa] */ - Scalar pressure(int phaseIdx) const + const Scalar& pressure(int phaseIdx) const { return pressure_[phaseIdx]; } /*! * \brief Set the pressure of a phase [Pa] */ - void setPressure(int phaseIdx, Scalar value) + void setPressure(int phaseIdx, const Scalar& value) { pressure_[phaseIdx] = value; } /*! @@ -68,8 +68,10 @@ public: template void assign(const FluidState& fs) { + typedef typename FluidState::Scalar FsScalar; + typedef Opm::MathToolbox FsToolbox; for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - pressure_[phaseIdx] = fs.pressure(phaseIdx); + pressure_[phaseIdx] = FsToolbox::template toLhs(fs.pressure(phaseIdx)); } } @@ -106,7 +108,7 @@ public: /*! * \brief The pressure of a fluid phase [Pa] */ - Scalar pressure(int phaseIdx) const + const Scalar& pressure(int phaseIdx) const { OPM_THROW(std::logic_error, "Pressure is not provided by this fluid state"); } diff --git a/opm/material/fluidstates/FluidStateSaturationModules.hpp b/opm/material/fluidstates/FluidStateSaturationModules.hpp index 148ee9895..f3b11851f 100644 --- a/opm/material/fluidstates/FluidStateSaturationModules.hpp +++ b/opm/material/fluidstates/FluidStateSaturationModules.hpp @@ -24,11 +24,12 @@ #ifndef OPM_FLUID_STATE_SATURATION_MODULES_HPP #define OPM_FLUID_STATE_SATURATION_MODULES_HPP -#include - #include #include +#include +#include + #include namespace Opm { @@ -51,13 +52,13 @@ public: /*! * \brief The saturation of a fluid phase [-] */ - Scalar saturation(int phaseIdx) const + const Scalar& saturation(int phaseIdx) const { return saturation_[phaseIdx]; } /*! * \brief Set the saturation of a phase [-] */ - void setSaturation(int phaseIdx, Scalar value) + void setSaturation(int phaseIdx, const Scalar& value) { saturation_[phaseIdx] = value; } /*! @@ -67,8 +68,10 @@ public: template void assign(const FluidState& fs) { + typedef typename FluidState::Scalar FsScalar; + typedef Opm::MathToolbox FsToolbox; for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - saturation_[phaseIdx] = fs.saturation(phaseIdx); + saturation_[phaseIdx] = FsToolbox::template toLhs(fs.saturation(phaseIdx)); } } @@ -105,7 +108,7 @@ public: /*! * \brief The saturation of a fluid phase [-] */ - Scalar saturation(int phaseIdx) const + const Scalar& saturation(int phaseIdx) const { OPM_THROW(std::runtime_error, "Saturation is not provided by this fluid state"); } /*! diff --git a/opm/material/fluidstates/FluidStateTemperatureModules.hpp b/opm/material/fluidstates/FluidStateTemperatureModules.hpp index 683b92f31..694f6b270 100644 --- a/opm/material/fluidstates/FluidStateTemperatureModules.hpp +++ b/opm/material/fluidstates/FluidStateTemperatureModules.hpp @@ -27,6 +27,7 @@ #include +#include #include #include @@ -125,11 +126,15 @@ public: template void assign(const FluidState& fs) { - temperature_ = fs.temperature(/*phaseIdx=*/0); + typedef Opm::MathToolbox Toolbox; + typedef Opm::MathToolbox FsToolbox; + + temperature_ = FsToolbox::template toLhs(fs.temperature(/*phaseIdx=*/0)); #ifndef NDEBUG for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - assert(std::abs(fs.temperature(phaseIdx) - temperature_) < 1e-30); + assert(std::abs(FsToolbox::value(fs.temperature(phaseIdx)) + - Toolbox::value(temperature_)) < 1e-30); } #endif } diff --git a/opm/material/fluidstates/FluidStateViscosityModules.hpp b/opm/material/fluidstates/FluidStateViscosityModules.hpp index 683ceb505..2c76316bd 100644 --- a/opm/material/fluidstates/FluidStateViscosityModules.hpp +++ b/opm/material/fluidstates/FluidStateViscosityModules.hpp @@ -25,15 +25,15 @@ #ifndef OPM_FLUID_STATE_VISCOSITY_MODULES_HPP #define OPM_FLUID_STATE_VISCOSITY_MODULES_HPP -#include - #include #include +#include +#include + #include namespace Opm { - /*! * \brief Module for the modular fluid state which stores the * viscosities explicitly. @@ -68,8 +68,10 @@ public: template void assign(const FluidState& fs) { + typedef typename FluidState::Scalar FsScalar; + typedef Opm::MathToolbox FsToolbox; for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { - viscosity_[phaseIdx] = fs.viscosity(phaseIdx); + viscosity_[phaseIdx] = FsToolbox::template toLhs(fs.viscosity(phaseIdx)); } } diff --git a/opm/material/fluidstates/PressureOverlayFluidState.hpp b/opm/material/fluidstates/PressureOverlayFluidState.hpp index 85f43cebd..91adf8b7a 100644 --- a/opm/material/fluidstates/PressureOverlayFluidState.hpp +++ b/opm/material/fluidstates/PressureOverlayFluidState.hpp @@ -34,10 +34,12 @@ namespace Opm { * pressures and takes all other quantities from an other * fluid state. */ -template +template class PressureOverlayFluidState { public: + typedef typename FluidState::Scalar Scalar; + enum { numPhases = FluidState::numPhases }; enum { numComponents = FluidState::numComponents }; diff --git a/opm/material/fluidstates/SaturationOverlayFluidState.hpp b/opm/material/fluidstates/SaturationOverlayFluidState.hpp index dd923508a..95e03cab6 100644 --- a/opm/material/fluidstates/SaturationOverlayFluidState.hpp +++ b/opm/material/fluidstates/SaturationOverlayFluidState.hpp @@ -34,10 +34,12 @@ namespace Opm { * saturations and takes all other quantities from an other * fluid state. */ -template +template class SaturationOverlayFluidState { public: + typedef typename FluidState::Scalar Scalar; + enum { numPhases = FluidState::numPhases }; enum { numComponents = FluidState::numComponents }; diff --git a/opm/material/fluidstates/TemperatureOverlayFluidState.hpp b/opm/material/fluidstates/TemperatureOverlayFluidState.hpp index 1ee9327dc..3c46d8624 100644 --- a/opm/material/fluidstates/TemperatureOverlayFluidState.hpp +++ b/opm/material/fluidstates/TemperatureOverlayFluidState.hpp @@ -32,10 +32,12 @@ namespace Opm { * temperatures and takes all other quantities from an other * fluid state. */ -template +template class TemperatureOverlayFluidState { public: + typedef typename FluidState::Scalar Scalar; + enum { numPhases = FluidState::numPhases }; enum { numComponents = FluidState::numComponents }; diff --git a/tests/checkFluidSystem.hpp b/tests/checkFluidSystem.hpp index af8f382f9..7697072ec 100644 --- a/tests/checkFluidSystem.hpp +++ b/tests/checkFluidSystem.hpp @@ -51,9 +51,9 @@ * \brief This is a fluid state which makes sure that only the quantities * allowed are accessed. */ -template > + class BaseFluidState = Opm::CompositionalFluidState > class HairSplittingFluidState : protected BaseFluidState { @@ -61,6 +61,8 @@ public: enum { numPhases = FluidSystem::numPhases }; enum { numComponents = FluidSystem::numComponents }; + typedef ScalarT Scalar; + HairSplittingFluidState() { // initially, do not allow anything @@ -221,6 +223,11 @@ void checkFluidState(const BaseFluidState &fs) // a fluid state must provide a checkDefined() method fs.checkDefined(); + // fluid states must export the types which they use as Scalars + typedef typename BaseFluidState::Scalar FsScalar; + static_assert(std::is_same::value, + "Fluid states must export the type they are given as scalar in an unmodified way"); + // make sure the fluid state provides all mandatory methods while (false) { Scalar OPM_UNUSED val; diff --git a/tests/test_fluidsystems.cpp b/tests/test_fluidsystems.cpp index 58622a4fb..ae381de5b 100644 --- a/tests/test_fluidsystems.cpp +++ b/tests/test_fluidsystems.cpp @@ -108,15 +108,15 @@ int main(int argc, char **argv) BaseFluidState baseFs; // TemperatureOverlayFluidState - { Opm::TemperatureOverlayFluidState fs(baseFs); + { Opm::TemperatureOverlayFluidState fs(baseFs); checkFluidState(fs); } // PressureOverlayFluidState - { Opm::PressureOverlayFluidState fs(baseFs); + { Opm::PressureOverlayFluidState fs(baseFs); checkFluidState(fs); } // SaturationOverlayFluidState - { Opm::SaturationOverlayFluidState fs(baseFs); + { Opm::SaturationOverlayFluidState fs(baseFs); checkFluidState(fs); } }