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

@@ -65,7 +65,7 @@ public:
     static char *phaseName(int phaseIdx)
     {
         OPM_THROW(std::runtime_error,
-                   "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' does not provide a phaseName() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' 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<Implementation>() << "' does not provide a isLiquid() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' 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<Implementation>() << "' does not provide a isIdealMixture() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' 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<Implementation>() << "' does not provide a isCompressible() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' 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<Implementation>() << "' does not provide a isIdealGas() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' 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<Implementation>() << "' does not provide a componentName() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' 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<Implementation>() << "' does not provide a molarMass() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' does not provide a molarMass() method!");
     }
 
     /*!
@@ -160,13 +160,13 @@ public:
      * \copydoc Doxygen::fluidSystemBaseParams
      * \copydoc Doxygen::phaseIdxParam
      */
-    template <class FluidState, class ParameterCache>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParameterCache = NullParameterCache>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
         OPM_THROW(std::runtime_error,
-                   "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' does not provide a density() method!");
+                  "Not implemented: The fluid system '" << Opm::className<Implementation>() << "' does not provide a density() method!");
     }
 
     /*!
@@ -183,11 +183,11 @@ public:
      * \copydoc Doxygen::phaseIdxParam
      * \copydoc Doxygen::compIdxParam
      */
-    template <class FluidState, class ParameterCache>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParameterCache = NullParameterCache>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
         OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className<Implementation>() << "'  does not provide a fugacityCoefficient() method!");
     }
@@ -198,10 +198,10 @@ public:
      * \copydoc Doxygen::fluidSystemBaseParams
      * \copydoc Doxygen::phaseIdxParam
      */
-    template <class FluidState, class ParameterCache>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParameterCache = NullParameterCache>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
         OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className<Implementation>() << "'  does not provide a viscosity() method!");
     }
@@ -223,13 +223,13 @@ public:
      * \copydoc Doxygen::phaseIdxParam
      * \copydoc Doxygen::compIdxParam
      */
-    template <class FluidState, class ParameterCache>
-    static Scalar diffusionCoefficient(const FluidState &fluidState,
-                                       const ParameterCache &paramCache,
-                                       int phaseIdx,
-                                       int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParameterCache = NullParameterCache>
+    static LhsEval diffusionCoefficient(const FluidState &fluidState,
+                                        const ParameterCache &paramCache,
+                                        int phaseIdx,
+                                        int compIdx)
     {
-         OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className<Implementation>() << "'  does not provide a diffusionCoefficient() method!");
+        OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className<Implementation>() << "'  does not provide a diffusionCoefficient() method!");
     }
 
     /*!
@@ -239,10 +239,10 @@ public:
      * \copydoc Doxygen::fluidSystemBaseParams
      * \copydoc Doxygen::phaseIdxParam
      */
-    template <class FluidState, class ParameterCache>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           const ParameterCache &paramCache,
-                           int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParameterCache = NullParameterCache>
+    static LhsEval enthalpy(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
         OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className<Implementation>() << "'  does not provide an enthalpy() method!");
     }
@@ -253,10 +253,10 @@ public:
      * \copydoc Doxygen::fluidSystemBaseParams
      * \copydoc Doxygen::phaseIdxParam
      */
-    template <class FluidState, class ParameterCache>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParameterCache = NullParameterCache>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx)
     {
         OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className<Implementation>() << "'  does not provide a thermalConductivity() method!");
     }
@@ -267,10 +267,10 @@ public:
      * \copydoc Doxygen::fluidSystemBaseParams
      * \copydoc Doxygen::phaseIdxParam
      */
-    template <class FluidState, class ParameterCache>
-    static Scalar heatCapacity(const FluidState &fluidState,
-                               const ParameterCache &paramCache,
-                               int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParameterCache = NullParameterCache>
+    static LhsEval heatCapacity(const FluidState &fluidState,
+                                const ParameterCache &paramCache,
+                                int phaseIdx)
     {
         OPM_THROW(std::runtime_error, "Not implemented: The fluid system '" << Opm::className<Implementation>() << "'  does not provide a heatCapacity() method!");
     }

opm/material/fluidsystems/BlackOilFluidSystem.hpp

@@ -38,27 +38,17 @@
 #include <array>
 
 namespace Opm {
-template <class Scalar>
-class OilPvtInterface;
-
-template <class Scalar>
-class GasPvtInterface;
-
-template <class Scalar>
-class WaterPvtInterface;
-
 namespace FluidSystems {
-
 /*!
  * \brief A fluid system which uses the black-oil parameters
  *        to calculate termodynamically meaningful quantities.
  */
-template <class Scalar>
-class BlackOil : public BaseFluidSystem<Scalar, BlackOil<Scalar> >
+template <class Scalar, class Evaluation = Scalar>
+class BlackOil : public BaseFluidSystem<Scalar, BlackOil<Scalar, Evaluation> >
 {
-    typedef Opm::GasPvtInterface<Scalar> GasPvtInterface;
-    typedef Opm::OilPvtInterface<Scalar> OilPvtInterface;
-    typedef Opm::WaterPvtInterface<Scalar> WaterPvtInterface;
+    typedef Opm::GasPvtInterface<Scalar, Evaluation> GasPvtInterface;
+    typedef Opm::OilPvtInterface<Scalar, Evaluation> OilPvtInterface;
+    typedef Opm::WaterPvtInterface<Scalar, Evaluation> WaterPvtInterface;
 
 public:
     //! \copydoc BaseFluidSystem::ParameterCache
@@ -286,26 +276,29 @@ public:
      * thermodynamic relations
      ****************************************/
     //! \copydoc BaseFluidSystem::density
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          ParameterCache &paramCache,
-                          const int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           ParameterCache &paramCache,
+                           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<FsEval> FsToolbox;
+
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
         int regionIdx = paramCache.regionIndex();
 
         switch (phaseIdx) {
-        case waterPhaseIdx: return waterDensity(T, p, regionIdx);
+        case waterPhaseIdx: return waterDensity<LhsEval>(T, p, regionIdx);
         case gasPhaseIdx: {
-            Scalar XgO = fluidState.massFraction(gasPhaseIdx, oilCompIdx);
-            return gasDensity(T, p, XgO, regionIdx);
+            const auto& XgO = FsToolbox::template toLhs<LhsEval>(fluidState.massFraction(gasPhaseIdx, oilCompIdx));
+            return gasDensity<LhsEval>(T, p, XgO, regionIdx);
         }
         case oilPhaseIdx: {
-            Scalar XoG = fluidState.massFraction(oilPhaseIdx, gasCompIdx);
-            return oilDensity(T, p, XoG, regionIdx);
+            const auto& XoG = FsToolbox::template toLhs<LhsEval>(fluidState.massFraction(oilPhaseIdx, gasCompIdx));
+            return oilDensity<LhsEval>(T, p, XoG, regionIdx);
         }
         }
 
@@ -313,49 +306,53 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::fugacityCoefficient
-    template <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       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<typename FluidState::Scalar> FsToolbox;
+
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        const auto& T = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(compIdx, T, p, regionIdx);
+        case gasPhaseIdx: return fugCoefficientInGas<LhsEval>(compIdx, T, p, regionIdx);
+        case oilPhaseIdx: return fugCoefficientInOil<LhsEval>(compIdx, T, p, regionIdx);
         }
 
         OPM_THROW(std::logic_error, "Unhandled phase or component index");
     }
 
     //! \copydoc BaseFluidSystem::viscosity
-    template <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx <= numPhases);
 
-        Scalar p = fluidState.pressure(phaseIdx);
-        Scalar T = fluidState.temperature(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
         int regionIdx = paramCache.regionIndex();
 
         switch (phaseIdx) {
         case oilPhaseIdx: {
-            Scalar XoG = fluidState.massFraction(oilPhaseIdx, gasCompIdx);
+            const auto& XoG = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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<const Opm::GasPvtInterface<Scalar> > gasPvt_;
-    static std::shared_ptr<const Opm::OilPvtInterface<Scalar> > oilPvt_;
-    static std::shared_ptr<const Opm::WaterPvtInterface<Scalar> > waterPvt_;
+    static std::shared_ptr<const Opm::GasPvtInterface<Scalar, Evaluation> > gasPvt_;
+    static std::shared_ptr<const Opm::OilPvtInterface<Scalar, Evaluation> > oilPvt_;
+    static std::shared_ptr<const Opm::WaterPvtInterface<Scalar, Evaluation> > waterPvt_;
 
     static bool enableDissolvedGas_;
     static bool enableVaporizedOil_;
@@ -554,39 +607,39 @@ private:
     static std::vector<std::array<Scalar, /*numComponents=*/3> > molarMass_;
 };
 
-template <class Scalar>
+template <class Scalar, class Evaluation>
 const Scalar
-BlackOil<Scalar>::surfaceTemperature = 273.15 + 15.56; // [K]
+BlackOil<Scalar, Evaluation>::surfaceTemperature = 273.15 + 15.56; // [K]
 
-template <class Scalar>
+template <class Scalar, class Evaluation>
 const Scalar
-BlackOil<Scalar>::surfacePressure = 101325.0; // [Pa]
+BlackOil<Scalar, Evaluation>::surfacePressure = 101325.0; // [Pa]
 
-template <class Scalar>
-bool BlackOil<Scalar>::enableDissolvedGas_;
+template <class Scalar, class Evaluation>
+bool BlackOil<Scalar, Evaluation>::enableDissolvedGas_;
 
-template <class Scalar>
-bool BlackOil<Scalar>::enableVaporizedOil_;
+template <class Scalar, class Evaluation>
+bool BlackOil<Scalar, Evaluation>::enableVaporizedOil_;
 
-template <class Scalar>
-std::shared_ptr<const OilPvtInterface<Scalar> >
-BlackOil<Scalar>::oilPvt_;
+template <class Scalar, class Evaluation>
+std::shared_ptr<const OilPvtInterface<Scalar, Evaluation> >
+BlackOil<Scalar, Evaluation>::oilPvt_;
 
-template <class Scalar>
-std::shared_ptr<const GasPvtInterface<Scalar> >
-BlackOil<Scalar>::gasPvt_;
+template <class Scalar, class Evaluation>
+std::shared_ptr<const GasPvtInterface<Scalar, Evaluation> >
+BlackOil<Scalar, Evaluation>::gasPvt_;
 
-template <class Scalar>
-std::shared_ptr<const WaterPvtInterface<Scalar> >
-BlackOil<Scalar>::waterPvt_;
+template <class Scalar, class Evaluation>
+std::shared_ptr<const WaterPvtInterface<Scalar, Evaluation> >
+BlackOil<Scalar, Evaluation>::waterPvt_;
 
-template <class Scalar>
+template <class Scalar, class Evaluation>
 std::vector<std::array<Scalar, 3> >
-BlackOil<Scalar>::referenceDensity_;
+BlackOil<Scalar, Evaluation>::referenceDensity_;
 
-template <class Scalar>
+template <class Scalar, class Evaluation>
 std::vector<std::array<Scalar, 3> >
-BlackOil<Scalar>::molarMass_;
+BlackOil<Scalar, Evaluation>::molarMass_;
 }} // namespace Opm, FluidSystems
 
 #endif

opm/material/fluidsystems/BrineCO2FluidSystem.hpp

@@ -226,30 +226,33 @@ public:
     /*!
      * \copydoc BaseFluidSystem::density
      */
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+        typedef MathToolbox<LhsEval> LhsToolbox;
+
         assert(0 <= phaseIdx && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const LhsEval& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& pressure = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, BrineIdx))));
+            LhsEval xlCO2 = LhsToolbox::min(1.0, LhsToolbox::max(0.0,  FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(gasPhaseIdx, BrineIdx))));
+        LhsEval xgCO2 = LhsToolbox::min(1.0, LhsToolbox::max(0.0,  FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const LhsEval& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& pressure = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+        typedef MathToolbox<LhsEval> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& pressure = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar diffusionCoefficient(const FluidState &fluidState,
-                                       const ParameterCache &paramCache,
-                                       int phaseIdx,
-                                       int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval diffusionCoefficient(const FluidState &fluidState,
+                                        const ParameterCache &paramCache,
+                                        int phaseIdx,
+                                        int compIdx)
     {
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const LhsEval& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& pressure = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
         if (phaseIdx == liquidPhaseIdx)
             return BinaryCoeffBrineCO2::liquidDiffCoeff(temperature, pressure);
 
@@ -376,30 +386,33 @@ public:
     /*!
      * \copydoc BaseFluidSystem::enthalpy
      */
-    template <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           const ParameterCache &paramCache,
-                           int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+        typedef MathToolbox<LhsEval> LhsToolbox;
+
         assert(0 <= phaseIdx && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const LhsEval& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& pressure = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.massFraction(gasPhaseIdx, CO2Idx));
+            const LhsEval& XBrine = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx)
     {
+        typedef MathToolbox<LhsEval> 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 <class FluidState>
-    static Scalar heatCapacity(const FluidState &fluidState,
-                               const ParameterCache &paramCache,
-                               int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval heatCapacity(const FluidState &fluidState,
+                                const ParameterCache &paramCache,
+                                int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        assert(0 <= phaseIdx && phaseIdx < numPhases);
+
+        const LhsEval& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& pressure = FsToolbox::template toLhs<LhsEval>(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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    static LhsEval liquidEnthalpyBrineCO2_(const LhsEval& T,
+                                           const LhsEval& p,
+                                           Scalar S, // salinity
+                                           const LhsEval& X_CO2_w)
     {
+        typedef MathToolbox<LhsEval> 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<Scalar>(i)) * std::pow(m, j);
             }
         }
         /* heat of dissolution for halite according to Michaelides 1971 */

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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    static Evaluation heatCapacity(const Evaluation& temperature, const Evaluation& pressure)
     { return Component::gasHeatCapacity(temperature, pressure); }
 };
 } // namespace Opm

opm/material/fluidsystems/H2OAirFluidSystem.hpp

@@ -249,17 +249,20 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::density
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+        typedef Opm::MathToolbox<LhsEval> LhsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p;
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        LhsEval p;
         if (isCompressible(phaseIdx))
-            p = fluidState.pressure(phaseIdx);
+            p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
+                const auto& xlAir = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.averageMolarMass(gasPhaseIdx))
+                    / LhsToolbox::max(1e-5, sumMoleFrac);
 
-            Scalar partialPressureH2O =
-                fluidState.moleFraction(gasPhaseIdx, H2OIdx)  *
-                fluidState.pressure(gasPhaseIdx);
+            LhsEval partialPressureH2O =
+                FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))
+                *FsToolbox::template toLhs<LhsEval>(fluidState.pressure(gasPhaseIdx));
 
-            Scalar partialPressureAir =
-                fluidState.moleFraction(gasPhaseIdx, AirIdx)  *
-                fluidState.pressure(gasPhaseIdx);
+            LhsEval partialPressureAir =
+                FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, AirIdx))
+                *FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
                             const ParameterCache &paramCache,
                             int phaseIdx)
     {
+        typedef Opm::MathToolbox<LhsEval> LhsToolbox;
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(phaseIdx, j))*phiIJ;
                     }
-                    muResult += fluidState.moleFraction(phaseIdx, i)*mu[i] / divisor;
+                    const auto& xAlphaI = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, i));
+                    muResult += xAlphaI*mu[i]/divisor;
                 }
                 return muResult;
             }
@@ -382,17 +386,19 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::fugacityCoefficient
-    template <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
                                       const ParameterCache &paramCache,
                                       int phaseIdx,
                                       int compIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
 
         if (phaseIdx == liquidPhaseIdx) {
             if (compIdx == H2OIdx)
@@ -406,14 +412,16 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::diffusionCoefficient
-    template <class FluidState>
-    static Scalar binaryDiffusionCoefficient(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval binaryDiffusionCoefficient(const FluidState &fluidState,
                                              const ParameterCache &paramCache,
                                              int phaseIdx,
                                              int compIdx)
     {
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
 
         if (phaseIdx == liquidPhaseIdx)
             return BinaryCoeff::H2O_Air::liquidDiffCoeff(T, p);
@@ -423,13 +431,15 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::enthalpy
-    template <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
                            const ParameterCache &paramCache,
                            int phaseIdx)
     {
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
 
             result +=
                 Air::gasEnthalpy(T, p) *
-                fluidState.massFraction(gasPhaseIdx, AirIdx);
+                FsToolbox::template toLhs<LhsEval>(fluidState.massFraction(gasPhaseIdx, AirIdx));
             return result;
         }
         OPM_THROW(std::logic_error, "Invalid phase index " << phaseIdx);
     }
 
     //! \copydoc BaseFluidSystem::thermalConductivity
-    template <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                          const ParameterCache &paramCache,
+                                          int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar temperature  = fluidState.temperature(phaseIdx) ;
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const LhsEval& temperature =
+            FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& pressure =
+            FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(phaseIdx, AirIdx));
+                const LhsEval& xH2O =
+                    FsToolbox::template toLhs<LhsEval>(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

opm/material/fluidsystems/H2OAirMesityleneFluidSystem.hpp

@@ -194,44 +194,46 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::density
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
-        Scalar T = fluidState.temperature(phaseIdx) ;
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const LhsEval& T = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx)) +
+                       Air::molarMass()*FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx)) +
+                       NAPL::molarMass()*FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.pressure(gasPhaseIdx));
+        const LhsEval& pH2O = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*pg;
+        const LhsEval& pAir = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*pg;
+        const LhsEval& pNAPL = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, NAPLIdx))*pg;
         return
             H2O::gasDensity(T, pH2O) +
             Air::gasDensity(T, pAir) +
@@ -239,13 +241,15 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::viscosity
-    template <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const LhsEval& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
+        const LhsEval& xgH2O = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
+        const LhsEval& xgNapl = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar diffusionCoefficient(const FluidState &fluidState,
-                                       const ParameterCache &paramCache,
-                                       int phaseIdx,
-                                       int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval diffusionCoefficient(const FluidState &fluidState,
+                                        const ParameterCache &paramCache,
+                                        int phaseIdx,
+                                        int compIdx)
     {
         return 0;
 #if 0
-        Scalar T = fluidState.temperature(phaseIdx) ;
-        Scalar p = fluidState.pressure(phaseIdx);
-        Scalar diffCont;
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const LhsEval& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
+            const LhsEval& xgw = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
+            const LhsEval& xgc = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
+            const LhsEval& xww = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
+            const LhsEval& xwc = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& p = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           const ParameterCache &paramCache,
-                           int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
-        Scalar T = fluidState.temperature(phaseIdx) ;
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const LhsEval& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const LhsEval& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
+            result += NAPL::gasEnthalpy(T, p) * FsToolbox::template toLhs<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
+            result += Air::gasEnthalpy(T, p) * FsToolbox::template toLhs<LhsEval>(fluidState.massFraction(gasPhaseIdx, NAPLIdx));
 
             return result;
         }
@@ -442,22 +444,26 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::thermalConductivity
-    template <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+            const LhsEval& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.temperature(phaseIdx));
+            const LhsEval& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+
+            return Air::gasThermalConductivity(T, p);
         }
 
         assert(phaseIdx == naplPhaseIdx);

opm/material/fluidsystems/H2OAirXyleneFluidSystem.hpp

@@ -161,62 +161,66 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::density
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
                           const ParameterCache &paramCache,
                           int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         if (phaseIdx == waterPhaseIdx) {
+            const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+            const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
+            const auto& xwAir = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
+            const auto& xwNapl = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+
+        const LhsEval& pH2O = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx))*p;
+        const LhsEval& pAir = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx))*p;
+        const LhsEval& pNAPL = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
                             const ParameterCache &paramCache,
                             int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
+        const auto& xgH2O = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
+        const auto& xgNapl = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar diffusionCoefficient(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval diffusionCoefficient(const FluidState &fluidState,
                                        const ParameterCache &paramCache,
                                        int phaseIdx,
                                        int compIdx)
     {
-        Scalar diffCont;
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+            const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(gasPhaseIdx, airIdx));
+            const LhsEval& xgw = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
+            const LhsEval& xgc = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(waterPhaseIdx, airIdx));
+            const LhsEval& xww = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(waterPhaseIdx, H2OIdx));
+            const LhsEval& xwc = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
                                       const ParameterCache &paramCache,
                                       int phaseIdx,
                                       int compIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
                            const ParameterCache &paramCache,
                            int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
+            result += hga * FsToolbox::template toLhs<LhsEval>(fluidState.massFraction(gasPhaseIdx, airIdx));
+            result += hgc * FsToolbox::template toLhs<LhsEval>(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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    static LhsEval NAPLPhaseDensity_(const LhsEval& T, const LhsEval& pn)
+    { return NAPL::liquidDensity(T, pn); }
 };
 } // namespace FluidSystems
 } // namespace Opm

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 <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<LhsEval> LhsToolbox;
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
 
-        Scalar sumMoleFrac = 0;
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
+                const auto& xlN2 = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.moleFraction(gasPhaseIdx, H2OIdx));
+        const auto& xgN2 = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<LhsEval> LhsToolbox;
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.moleFraction(phaseIdx, j))
+                        /sumx*phiIJ;
                 }
-                muResult += fluidState.moleFraction(phaseIdx, i)/sumx * mu[i] / divisor;
+                muResult +=
+                    FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, i))
+                    /sumx*mu[i]/divisor;
             }
             return muResult;
         }
     }
 
     //! \copydoc BaseFluidSystem::fugacityCoefficient
-    template <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       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<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar diffusionCoefficient(const FluidState &fluidState,
-                                             const ParameterCache &paramCache,
-                                             int phaseIdx,
-                                             int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval diffusionCoefficient(const FluidState &fluidState,
+                                        const ParameterCache &paramCache,
+                                        int phaseIdx,
+                                        int compIdx)
 
     {
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
 
         // liquid phase
         if (phaseIdx == liquidPhaseIdx)
@@ -416,13 +429,15 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::enthalpy
-    template <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           const ParameterCache &paramCache,
-                           int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.massFraction(gasPhaseIdx, H2OIdx));
+        const auto& XgN2 = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx) ;
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
+            const auto& xN2 = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar heatCapacity(const FluidState &fluidState,
-                               const ParameterCache &paramCache,
-                               int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval heatCapacity(const FluidState &fluidState,
+                                const ParameterCache &paramCache,
+                                int phaseIdx)
     {
-        if (phaseIdx == liquidPhaseIdx) {
-            return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx),
-                                           fluidState.pressure(phaseIdx));
-        }
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        const auto& xAlphaH2O = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, H2OIdx));
+        const auto& xAlphaN2 = FsToolbox::template toLhs<LhsEval>(fluidState.moleFraction(phaseIdx, N2Idx));
+        const auto& XAlphaH2O = FsToolbox::template toLhs<LhsEval>(fluidState.massFraction(phaseIdx, H2OIdx));
+        const auto& XAlphaN2 = FsToolbox::template toLhs<LhsEval>(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;
     }
 };
 

opm/material/fluidsystems/H2ON2LiquidPhaseFluidSystem.hpp

@@ -245,19 +245,21 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::density
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(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<LhsEval>(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<LhsEval>(fluidState.moleFraction(liquidPhaseIdx, H2OIdx));
+            const auto& xlN2 = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(phaseIdx == liquidPhaseIdx);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
 
         // assume pure water for the liquid phase
         return H2O::liquidViscosity(T, p);
     }
 
     //! \copydoc BaseFluidSystem::fugacityCoefficient
-    template <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
 
         if (compIdx == H2OIdx)
             return H2O::vaporPressure(T)/p;
@@ -315,31 +319,35 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::diffusionCoefficient
-    template <class FluidState>
-    static Scalar diffusionCoefficient(const FluidState &fluidState,
-                                       const ParameterCache &paramCache,
-                                       int phaseIdx,
-                                       int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval diffusionCoefficient(const FluidState &fluidState,
+                                        const ParameterCache &paramCache,
+                                        int phaseIdx,
+                                        int compIdx)
 
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(phaseIdx == liquidPhaseIdx);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
 
         return BinaryCoeff::H2O_N2::liquidDiffCoeff(T, p);
     }
 
     //! \copydoc BaseFluidSystem::enthalpy
-    template <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           const ParameterCache &paramCache,
-                           int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert (phaseIdx == liquidPhaseIdx);
 
-        Scalar T = fluidState.temperature(phaseIdx);
-        Scalar p = fluidState.pressure(phaseIdx);
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
         Valgrind::CheckDefined(T);
         Valgrind::CheckDefined(p);
 
@@ -348,32 +356,38 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::thermalConductivity
-    template <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      const int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       const int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+            const auto& p = FsToolbox::template toLhs<LhsEval>(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 <class FluidState>
-    static Scalar heatCapacity(const FluidState &fluidState,
-                               const ParameterCache &paramCache,
-                               int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval heatCapacity(const FluidState &fluidState,
+                                const ParameterCache &paramCache,
+                                int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert (phaseIdx == liquidPhaseIdx);
 
-        return H2O::liquidHeatCapacity(fluidState.temperature(phaseIdx),
-                                       fluidState.pressure(phaseIdx));
+        const auto& T = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+
+        return H2O::liquidHeatCapacity(T, p);
     }
 };
 

opm/material/fluidsystems/LiquidPhase.hpp

@@ -74,35 +74,43 @@ public:
     { return Component::triplePressure(); }
 
     //! \copydoc GasPhase::vaporPressure
-    static Scalar vaporPressure(Scalar temperature)
+    template <class Evaluation>
+    static Evaluation vaporPressure(const Evaluation& temperature)
     {  return Component::vaporPressure(temperature); }
 
     //! \copydoc GasPhase::density
-    static Scalar density(Scalar temperature, Scalar pressure)
+    template <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    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 <class Evaluation>
+    static Evaluation heatCapacity(const Evaluation& temperature, const Evaluation& pressure)
     { return Component::liquidHeatCapacity(temperature, pressure); }
 };
 } // namespace Opm

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 <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        return Fluid::density(T, p);
     }
 
     //! \copydoc BaseFluidSystem::viscosity
-    template <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        return Fluid::viscosity(T, p);
     }
 
     //! \copydoc BaseFluidSystem::fugacityCoefficient
-    template <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
         assert(0 <= compIdx  && compIdx < numComponents);
@@ -223,42 +227,48 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::enthalpy
-    template <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           const ParameterCache &paramCache,
-                           int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        return Fluid::enthalpy(T, p);
     }
 
     //! \copydoc BaseFluidSystem::thermalConductivity
-    template <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        return Fluid::thermalConductivity(T, p);
     }
 
     //! \copydoc BaseFluidSystem::heatCapacity
-    template <class FluidState>
-    static Scalar heatCapacity(const FluidState &fluidState,
-                               const ParameterCache &paramCache,
-                               int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval heatCapacity(const FluidState &fluidState,
+                                const ParameterCache &paramCache,
+                                int phaseIdx)
     {
+        typedef Opm::MathToolbox<typename FluidState::Scalar> 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<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& p = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
+        return Fluid::heatCapacity(T, p);
     }
 };
 

opm/material/fluidsystems/Spe5FluidSystem.hpp

@@ -352,23 +352,27 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::density
-    template <class FluidState>
+    template <class FluidState, class Evaluation = Scalar>
     static Scalar density(const FluidState &fluidState,
                           const ParameterCache &paramCache,
                           int phaseIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
+        static_assert(std::is_same<Evaluation, Scalar>::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 <class FluidState>
+    template <class FluidState, class Evaluation = Scalar>
     static Scalar viscosity(const FluidState &fluidState,
                             const ParameterCache &paramCache,
                             int phaseIdx)
     {
         assert(0 <= phaseIdx  && phaseIdx <= numPhases);
+        static_assert(std::is_same<Evaluation, Scalar>::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 <class FluidState>
+    template <class FluidState, class Evaluation = Scalar>
     static Scalar fugacityCoefficient(const FluidState &fluidState,
                                       const ParameterCache &paramCache,
                                       int phaseIdx,
@@ -395,6 +399,8 @@ public:
     {
         assert(0 <= phaseIdx  && phaseIdx <= numPhases);
         assert(0 <= compIdx  && compIdx <= numComponents);
+        static_assert(std::is_same<Evaluation, Scalar>::value,
+                      "The SPE-5 fluid system is currently only implemented for the scalar case.");
 
         if (phaseIdx == oilPhaseIdx || phaseIdx == gasPhaseIdx)
             return PengRobinsonMixture::computeFugacityCoefficient(fluidState,

opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp

@@ -52,7 +52,7 @@ namespace FluidSystems {
  */
 template <class Scalar, class WettingPhase, class NonwettingPhase>
 class TwoPhaseImmiscible
-: public BaseFluidSystem<Scalar, TwoPhaseImmiscible<Scalar, WettingPhase, NonwettingPhase> >
+    : public BaseFluidSystem<Scalar, TwoPhaseImmiscible<Scalar, WettingPhase, NonwettingPhase> >
 {
     // do not try to instanciate this class, it has only static members!
     TwoPhaseImmiscible()
@@ -216,42 +216,48 @@ public:
     }
 
     //! \copydoc BaseFluidSystem::density
-    template <class FluidState>
-    static Scalar density(const FluidState &fluidState,
-                          const ParameterCache &paramCache,
-                          int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval density(const FluidState &fluidState,
+                           const ParameterCache &paramCache,
+                           int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const auto& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& pressure = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
         if (phaseIdx == wettingPhaseIdx)
             return WettingPhase::density(temperature, pressure);
         return NonwettingPhase::density(temperature, pressure);
     }
 
     //! \copydoc BaseFluidSystem::viscosity
-    template <class FluidState>
-    static Scalar viscosity(const FluidState &fluidState,
-                            const ParameterCache &paramCache,
-                            int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval viscosity(const FluidState &fluidState,
+                             const ParameterCache &paramCache,
+                             int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const auto& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& pressure = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
         if (phaseIdx == wettingPhaseIdx)
             return WettingPhase::viscosity(temperature, pressure);
         return NonwettingPhase::viscosity(temperature, pressure);
     }
 
     //! \copydoc BaseFluidSystem::fugacityCoefficient
-    template <class FluidState>
-    static Scalar fugacityCoefficient(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx,
-                                      int compIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval fugacityCoefficient(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx,
+                                       int compIdx)
     {
+        typedef MathToolbox<LhsEval> 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<Scalar>::infinity();
+            return LhsToolbox::createConstant(1.0);
+        return LhsToolbox::createConstant(std::numeric_limits<Scalar>::infinity());
     }
 
     //! \copydoc BaseFluidSystem::enthalpy
-    template <class FluidState>
-    static Scalar enthalpy(const FluidState &fluidState,
-                           const ParameterCache &paramCache,
-                           int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval enthalpy(const FluidState &fluidState,
+                            const ParameterCache &paramCache,
+                            int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const auto& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& pressure = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
         if (phaseIdx == wettingPhaseIdx)
             return WettingPhase::enthalpy(temperature, pressure);
         return NonwettingPhase::enthalpy(temperature, pressure);
     }
 
     //! \copydoc BaseFluidSystem::thermalConductivity
-    template <class FluidState>
-    static Scalar thermalConductivity(const FluidState &fluidState,
-                                      const ParameterCache &paramCache,
-                                      int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval thermalConductivity(const FluidState &fluidState,
+                                       const ParameterCache &paramCache,
+                                       int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const auto& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& pressure = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
         if (phaseIdx == wettingPhaseIdx)
             return WettingPhase::thermalConductivity(temperature, pressure);
         return NonwettingPhase::thermalConductivity(temperature, pressure);
     }
 
     //! \copydoc BaseFluidSystem::heatCapacity
-    template <class FluidState>
-    static Scalar heatCapacity(const FluidState &fluidState,
-                               const ParameterCache &paramCache,
-                               int phaseIdx)
+    template <class FluidState, class LhsEval = typename FluidState::Scalar>
+    static LhsEval heatCapacity(const FluidState &fluidState,
+                                const ParameterCache &paramCache,
+                                int phaseIdx)
     {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
         assert(0 <= phaseIdx  && phaseIdx < numPhases);
 
-        Scalar temperature = fluidState.temperature(phaseIdx);
-        Scalar pressure = fluidState.pressure(phaseIdx);
+        const auto& temperature = FsToolbox::template toLhs<LhsEval>(fluidState.temperature(phaseIdx));
+        const auto& pressure = FsToolbox::template toLhs<LhsEval>(fluidState.pressure(phaseIdx));
         if (phaseIdx == wettingPhaseIdx)
             return WettingPhase::heatCapacity(temperature, pressure);
         return NonwettingPhase::heatCapacity(temperature, pressure);

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 Scalar>
-class ConstantCompressibilityOilPvt : public OilPvtInterface<Scalar>
+template <class Scalar, class Evaluation = Scalar>
+class ConstantCompressibilityOilPvt  : public OilPvtInterfaceTemplateWrapper<Scalar,
+                                                                             Evaluation,
+                                                                             ConstantCompressibilityOilPvt<Scalar, Evaluation> >
 {
-    typedef FluidSystems::BlackOil<Scalar> BlackOilFluidSystem;
+    friend class OilPvtInterfaceTemplateWrapper<Scalar,
+                                                Evaluation,
+                                                ConstantCompressibilityOilPvt<Scalar, Evaluation> >;
+
+    typedef FluidSystems::BlackOil<Scalar, Evaluation> BlackOilFluidSystem;
 
     typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
     typedef std::vector<std::pair<Scalar, Scalar> > 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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    LhsEval gasDissolutionFactor_(int regionIdx,
+                                  const LhsEval& temperature,
+                                  const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval oilSaturationPressure_(int regionIdx,
+                                   const LhsEval& temperature,
+                                   const LhsEval& XoG) const
+    {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval saturatedOilGasMassFraction_(int regionIdx,
+                                         const LhsEval& temperature,
+                                         const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> Toolbox;
+
+        return Toolbox::createConstant(0.0); /* this is dead oil! */
+    }
+
+    template <class LhsEval>
+    LhsEval saturatedOilGasMoleFraction_(int regionIdx,
+                                         const LhsEval& temperature,
+                                         const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> Toolbox;
+
+        return Toolbox::createConstant(0.0); /* this is dead oil! */
+    }
 
 private:
     std::vector<Scalar> oilReferencePressure_;

opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp

@@ -26,7 +26,6 @@
 #include "WaterPvtInterface.hpp"
 
 #include <opm/material/common/OpmFinal.hpp>
-
 #if HAVE_OPM_PARSER
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #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 Scalar>
-class ConstantCompressibilityWaterPvt : public WaterPvtInterface<Scalar>
+template <class Scalar, class Evaluation = Scalar>
+class ConstantCompressibilityWaterPvt
+    : public WaterPvtInterfaceTemplateWrapper<Scalar, Evaluation, ConstantCompressibilityWaterPvt<Scalar, Evaluation> >
 {
-    typedef FluidSystems::BlackOil<Scalar> BlackOilFluidSystem;
+    friend class WaterPvtInterfaceTemplateWrapper<Scalar, Evaluation, ConstantCompressibilityWaterPvt<Scalar, Evaluation> >;
+
+    typedef FluidSystems::BlackOil<Scalar, Evaluation> BlackOilFluidSystem;
 
     typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
     typedef std::vector<std::pair<Scalar, Scalar> > 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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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

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 Scalar>
-class DeadOilPvt : public OilPvtInterface<Scalar>
+template <class Scalar, class Evaluation = Scalar>
+class DeadOilPvt : public OilPvtInterfaceTemplateWrapper<Scalar,
+                                                         Evaluation,
+                                                         DeadOilPvt<Scalar, Evaluation> >
 {
-    typedef FluidSystems::BlackOil<Scalar> BlackOilFluidSystem;
+    friend class OilPvtInterfaceTemplateWrapper<Scalar,
+                                                Evaluation,
+                                                DeadOilPvt<Scalar, Evaluation> >;
+
+    typedef FluidSystems::BlackOil<Scalar, Evaluation> BlackOilFluidSystem;
 
     typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
     typedef std::vector<std::pair<Scalar, Scalar> > SamplingPoints;
@@ -96,8 +102,8 @@ public:
     {
         assert(pvdoTable.numRows() > 1);
 
-        const std::vector<Scalar>& BColumn(pvdoTable.getFormationFactorColumn());
-        std::vector<Scalar> invBColumn(pvdoTable.getFormationFactorColumn());
+        const auto& BColumn(pvdoTable.getFormationFactorColumn());
+        std::vector<Scalar> 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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    LhsEval gasDissolutionFactor_(int regionIdx,
+                                  const LhsEval& temperature,
+                                  const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval oilSaturationPressure_(int regionIdx,
+                                   const LhsEval& temperature,
+                                   const LhsEval& XoG) const
+    {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval saturatedOilGasMassFraction_(int regionIdx,
+                                         const LhsEval& temperature,
+                                         const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> Toolbox;
+
+        return Toolbox::createConstant(0.0); /* this is dead oil! */
+    }
+
+    template <class LhsEval>
+    LhsEval saturatedOilGasMoleFraction_(int regionIdx,
+                                         const LhsEval& temperature,
+                                         const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> Toolbox;
+
+        return Toolbox::createConstant(0.0); /* this is dead oil! */
+    }
 
 private:
     std::vector<TabulatedOneDFunction> inverseOilB_;

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 Scalar>
-class DryGasPvt : public GasPvtInterface<Scalar>
+template <class Scalar, class Evaluation = Scalar>
+class DryGasPvt
+    : public GasPvtInterfaceTemplateWrapper<Scalar, Evaluation, DryGasPvt<Scalar, Evaluation> >
 {
-    typedef FluidSystems::BlackOil<Scalar> BlackOilFluidSystem;
+    friend class GasPvtInterfaceTemplateWrapper<Scalar, Evaluation, DryGasPvt<Scalar, Evaluation> >;
+
+    typedef FluidSystems::BlackOil<Scalar, Evaluation> BlackOilFluidSystem;
 
     typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
     typedef std::vector<std::pair<Scalar, Scalar> > 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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    LhsEval fugacityCoefficient_(int regionIdx,
+                                 const LhsEval& temperature,
+                                 const LhsEval& pressure,
+                                 int compIdx) const
     {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval gasSaturationPressure_(int regionIdx,
+                                   const LhsEval& temperature,
+                                   const LhsEval& XgO) const
+    {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval oilVaporizationFactor_(int regionIdx,
+                                   const LhsEval& temperature,
+                                   const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval saturatedGasOilMassFraction_(int regionIdx,
+                                         const LhsEval& temperature,
+                                         const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    LhsEval saturatedGasOilMoleFraction_(int regionIdx,
+                                         const LhsEval& temperature,
+                                         const LhsEval& pressure) const
+    {
+        typedef Opm::MathToolbox<LhsEval> Toolbox;
+        return Toolbox::createConstant(0.0);  // this is dry gas!
+    }
 
 private:
     std::vector<TabulatedOneDFunction> inverseGasB_;

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 <opm/material/common/OpmFinal.hpp>
 
+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 <class Scalar>
+template <class Scalar, class Evaluation = Scalar>
 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 Scalar>
+class GasPvtInterface<Scalar, Scalar>
+{
+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 Scalar, class Evaluation, class Implementation>
+class GasPvtInterfaceTemplateWrapper : public GasPvtInterface<Scalar, Evaluation>
+{
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(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 Scalar, class Implementation>
+class GasPvtInterfaceTemplateWrapper<Scalar, Scalar, Implementation>
+    : public GasPvtInterface<Scalar, Scalar>
+{
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(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

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 Scalar>
-class LiveOilPvt : public OilPvtInterface<Scalar>
+template <class Scalar, class Evaluation = Scalar>
+class LiveOilPvt : public OilPvtInterfaceTemplateWrapper<Scalar, Evaluation, LiveOilPvt<Scalar, Evaluation> >
 {
-    typedef FluidSystems::BlackOil<Scalar> BlackOilFluidSystem;
+    friend class OilPvtInterfaceTemplateWrapper<Scalar, Evaluation, LiveOilPvt<Scalar, Evaluation> >;
+
+    typedef FluidSystems::BlackOil<Scalar, Evaluation> BlackOilFluidSystem;
 
     typedef Opm::UniformXTabulated2DFunction<Scalar> TabulatedTwoDFunction;
     typedef Opm::Tabulated1DFunction<Scalar> TabulatedOneDFunction;
@@ -61,11 +63,11 @@ class LiveOilPvt : public OilPvtInterface<Scalar>
 public:
     void setNumRegions(int numRegions)
     {
-        if (static_cast<int>(inverseOilB_.size()) < numRegions) {
-            inverseOilB_.resize(numRegions);
-            inverseOilBMu_.resize(numRegions);
-            oilMu_.resize(numRegions);
-            gasDissolutionFactor_.resize(numRegions);
+        if (static_cast<int>(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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    LhsEval oilSaturationPressure_(int regionIdx,
+                                   const LhsEval& temperature,
+                                   const LhsEval& XoG) const
     {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    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 <class LhsEval>
+    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<Scalar, Scalar> val(XoG, pSat);
             pSatSamplePoints.push_back(val);
@@ -527,10 +546,10 @@ private:
                                                                   /*type=*/Spline::Monotonic);
     }
 
-    std::vector<TabulatedTwoDFunction> inverseOilB_;
-    std::vector<TabulatedTwoDFunction> oilMu_;
-    std::vector<TabulatedTwoDFunction> inverseOilBMu_;
-    std::vector<TabulatedOneDFunction> gasDissolutionFactor_;
+    std::vector<TabulatedTwoDFunction> inverseOilBTable_;
+    std::vector<TabulatedTwoDFunction> oilMuTable_;
+    std::vector<TabulatedTwoDFunction> inverseOilBMuTable_;
+    std::vector<TabulatedOneDFunction> gasDissolutionFactorTable_;
     std::vector<Spline> saturationPressureSpline_;
 };
 

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 <opm/material/common/OpmFinal.hpp>
 
+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 <class Scalar>
+template <class Scalar, class Evaluation = Scalar>
 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 Scalar>
+class OilPvtInterface<Scalar, Scalar>
+{
 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 Scalar, class Evaluation, class Implementation>
+class OilPvtInterfaceTemplateWrapper
+    : public OilPvtInterface<Scalar, Evaluation>
+{
+    Evaluation viscosity(int regionIdx,
+                         const Evaluation& temperature,
+                         const Evaluation& pressure,
+                         const Evaluation& XoG) const OPM_FINAL
+    { return asImp_().template viscosity_<Evaluation>(regionIdx, temperature, pressure, XoG); }
+    Scalar viscosity(int regionIdx,
+                     Scalar temperature,
+                     Scalar pressure,
+                     Scalar XoG) const OPM_FINAL
+    { return asImp_().template viscosity_<Scalar>(regionIdx, temperature, pressure, XoG); }
+
+    Evaluation formationVolumeFactor(int regionIdx,
+                                     const Evaluation& temperature,
+                                     const Evaluation& pressure,
+                                     const Evaluation& XoG) const OPM_FINAL
+    { return asImp_().template formationVolumeFactor_<Evaluation>(regionIdx, temperature, pressure, XoG); }
+    Scalar formationVolumeFactor(int regionIdx,
+                                 Scalar temperature,
+                                 Scalar pressure,
+                                 Scalar XoG) const OPM_FINAL
+    { return asImp_().template formationVolumeFactor_<Scalar>(regionIdx, temperature, pressure, XoG); }
+
+    Evaluation density(int regionIdx,
+                       const Evaluation& temperature,
+                       const Evaluation& pressure,
+                       const Evaluation& XoG) const OPM_FINAL
+    { return asImp_().template density_<Evaluation>(regionIdx, temperature, pressure, XoG); }
+    Scalar density(int regionIdx,
+                   Scalar temperature,
+                   Scalar pressure,
+                   Scalar XoG) const OPM_FINAL
+    { return asImp_().template density_<Scalar>(regionIdx, temperature, pressure, XoG); }
+
+    Evaluation fugacityCoefficient(int regionIdx,
+                                   const Evaluation& temperature,
+                                   const Evaluation& pressure,
+                                   int compIdx) const OPM_FINAL
+    { return asImp_().template fugacityCoefficient_<Evaluation>(regionIdx, temperature, pressure, compIdx); }
+    Scalar fugacityCoefficient(int regionIdx,
+                               Scalar temperature,
+                               Scalar pressure,
+                               int compIdx) const OPM_FINAL
+    { return asImp_().template fugacityCoefficient_<Scalar>(regionIdx, temperature, pressure, compIdx); }
+
+    Evaluation gasDissolutionFactor(int regionIdx,
+                                    const Evaluation& temperature,
+                                    const Evaluation& pressure) const OPM_FINAL
+    { return asImp_().template gasDissolutionFactor_<Evaluation>(regionIdx, temperature, pressure); }
+    Scalar gasDissolutionFactor(int regionIdx,
+                                Scalar temperature,
+                                Scalar pressure) const OPM_FINAL
+    { return asImp_().template gasDissolutionFactor_<Scalar>(regionIdx, temperature, pressure); }
+
+    Evaluation oilSaturationPressure(int regionIdx,
+                                     const Evaluation& temperature,
+                                     const Evaluation& XoG) const OPM_FINAL
+    { return asImp_().template oilSaturationPressure_<Evaluation>(regionIdx, temperature, XoG); }
+    Scalar oilSaturationPressure(int regionIdx,
+                                 Scalar temperature,
+                                 Scalar XoG) const OPM_FINAL
+    { return asImp_().template oilSaturationPressure_<Scalar>(regionIdx, temperature, XoG); }
+
+    Evaluation saturatedOilGasMassFraction(int regionIdx,
+                                           const Evaluation& temperature,
+                                           const Evaluation& pressure) const OPM_FINAL
+    { return asImp_().template saturatedOilGasMassFraction_<Evaluation>(regionIdx, temperature, pressure); }
+    Scalar saturatedOilGasMassFraction(int regionIdx,
+                                       Scalar temperature,
+                                       Scalar pressure) const OPM_FINAL
+    { return asImp_().template saturatedOilGasMassFraction_<Scalar>(regionIdx, temperature, pressure); }
+
+    Evaluation saturatedOilGasMoleFraction(int regionIdx,
+                                           const Evaluation& temperature,
+                                           const Evaluation& pressure) const OPM_FINAL
+    { return asImp_().template saturatedOilGasMoleFraction_<Evaluation>(regionIdx, temperature, pressure); }
+    Scalar saturatedOilGasMoleFraction(int regionIdx,
+                                       Scalar temperature,
+                                       Scalar pressure) const OPM_FINAL
+    { return asImp_().template saturatedOilGasMoleFraction_<Scalar>(regionIdx, temperature, pressure); }
+
+private:
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+
+template <class Scalar, class Implementation>
+class OilPvtInterfaceTemplateWrapper<Scalar, Scalar, Implementation>
+    : public OilPvtInterface<Scalar, Scalar>
+{
+    Scalar viscosity(int regionIdx,
+                     Scalar temperature,
+                     Scalar pressure,
+                     Scalar XoG) const OPM_FINAL
+    { return asImp_().template viscosity_<Scalar>(regionIdx, temperature, pressure, XoG); }
+    Scalar formationVolumeFactor(int regionIdx,
+                                 Scalar temperature,
+                                 Scalar pressure,
+                                 Scalar XoG) const OPM_FINAL
+    { return asImp_().template formationVolumeFactor_<Scalar>(regionIdx, temperature, pressure, XoG); }
+    Scalar density(int regionIdx,
+                   Scalar temperature,
+                   Scalar pressure,
+                   Scalar XoG) const OPM_FINAL
+    { return asImp_().template density_<Scalar>(regionIdx, temperature, pressure, XoG); }
+    Scalar fugacityCoefficient(int regionIdx,
+                               Scalar temperature,
+                               Scalar pressure,
+                               int compIdx) const OPM_FINAL
+    { return asImp_().template fugacityCoefficient_<Scalar>(regionIdx, temperature, pressure, compIdx); }
+    Scalar gasDissolutionFactor(int regionIdx,
+                                Scalar temperature,
+                                Scalar pressure) const OPM_FINAL
+    { return asImp_().template gasDissolutionFactor_<Scalar>(regionIdx, temperature, pressure); }
+    Scalar oilSaturationPressure(int regionIdx,
+                                 Scalar temperature,
+                                 Scalar XoG) const OPM_FINAL
+    { return asImp_().template oilSaturationPressure_<Scalar>(regionIdx, temperature, XoG); }
+    Scalar saturatedOilGasMassFraction(int regionIdx,
+                                       Scalar temperature,
+                                       Scalar pressure) const OPM_FINAL
+    { return asImp_().template saturatedOilGasMassFraction_<Scalar>(regionIdx, temperature, pressure); }
+    Scalar saturatedOilGasMoleFraction(int regionIdx,
+                                       Scalar temperature,
+                                       Scalar pressure) const OPM_FINAL
+    { return asImp_().template saturatedOilGasMoleFraction_<Scalar>(regionIdx, temperature, pressure); }
+
+private:
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+
+
 } // namespace Opm
 
 #endif

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 <opm/material/common/OpmFinal.hpp>
 
+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 <class Scalar>
+template <class Scalar, class Evaluation = Scalar>
 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 Scalar>
+class WaterPvtInterface<Scalar, Scalar>
+{
+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 Scalar, class Evaluation, class Implementation>
+class WaterPvtInterfaceTemplateWrapper : public WaterPvtInterface<Scalar, Evaluation>
+{
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(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 Scalar, class Implementation>
+class WaterPvtInterfaceTemplateWrapper<Scalar, Scalar, Implementation>
+    : public WaterPvtInterface<Scalar, Scalar>
+{
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(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

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 Scalar>
-class WetGasPvt : public GasPvtInterface<Scalar>
+template <class Scalar, class Evaluation = Scalar>
+class WetGasPvt
+    : public GasPvtInterfaceTemplateWrapper<Scalar, Evaluation, WetGasPvt<Scalar, Evaluation> >
 {
-    typedef FluidSystems::BlackOil<Scalar> BlackOilFluidSystem;
+    friend class GasPvtInterfaceTemplateWrapper<Scalar, Evaluation, WetGasPvt<Scalar, Evaluation> >;
+
+    typedef FluidSystems::BlackOil<Scalar, Evaluation> BlackOilFluidSystem;
 
     typedef Opm::UniformXTabulated2DFunction<Scalar> TabulatedTwoDFunction;
     typedef Opm::Tabulated1DFunction<Scalar> 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<int>(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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    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 <class LhsEval>
+    LhsEval fugacityCoefficient_(int regionIdx,
+                                 const LhsEval& temperature,
+                                 const LhsEval& pressure,
+                                 int compIdx) const
     {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    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 <class LhsEval>
+    LhsEval gasSaturationPressure_(int regionIdx,
+                                   const LhsEval& temperature,
+                                   const LhsEval& XgO) const
     {
+        typedef Opm::MathToolbox<LhsEval> 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 <class LhsEval>
+    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 <class LhsEval>
+    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<Scalar, Scalar> val(XgO, pSat);
             pSatSamplePoints.push_back(val);
@@ -525,7 +539,7 @@ private:
     std::vector<TabulatedTwoDFunction> inverseGasB_;
     std::vector<TabulatedTwoDFunction> gasMu_;
     std::vector<TabulatedTwoDFunction> inverseGasBMu_;
-    std::vector<TabulatedOneDFunction> oilVaporizationFactor_;
+    std::vector<TabulatedOneDFunction> oilVaporizationFactorTable_;
     std::vector<Spline> saturationPressureSpline_;
 };
 

tests/checkFluidSystem.hpp

@@ -253,7 +253,7 @@ void checkFluidState(const BaseFluidState &fs)
 /*!
  * \brief Checks whether a fluid system adheres to the specification.
  */
-template <class Scalar, class FluidSystem>
+template <class Scalar, class FluidSystem, class RhsEval, class LhsEval>
 void checkFluidSystem()
 {
     std::cout << "Testing fluid system '" << Opm::className<FluidSystem>() << "'\n";
@@ -263,7 +263,8 @@ void checkFluidSystem()
     enum { numPhases = FluidSystem::numPhases };
     enum { numComponents = FluidSystem::numComponents };
 
-    HairSplittingFluidState<Scalar, FluidSystem> fs;
+    typedef HairSplittingFluidState<RhsEval, FluidSystem> 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<decltype(tmpVal), RhsEval>::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {};
+        try { val = FluidSystem::template density<FluidState, LhsEval>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { scalarVal = FluidSystem::template density<FluidState, Scalar>(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<decltype(tmpVal), RhsEval>::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<decltype(tmpVal), RhsEval>::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<decltype(tmpVal), RhsEval>::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<decltype(tmpVal), RhsEval>::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {};
+        try { val = FluidSystem::template viscosity<FluidState, LhsEval>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { val = FluidSystem::template enthalpy<FluidState, LhsEval>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { val = FluidSystem::template heatCapacity<FluidState, LhsEval>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { val = FluidSystem::template thermalConductivity<FluidState, LhsEval>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { scalarVal = FluidSystem::template viscosity<FluidState, Scalar>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { scalarVal = FluidSystem::template enthalpy<FluidState, Scalar>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { scalarVal = FluidSystem::template heatCapacity<FluidState, Scalar>(fs, paramCache, phaseIdx); } catch (...) {};
+        try { scalarVal = FluidSystem::template thermalConductivity<FluidState, Scalar>(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<decltype(tmpVal), RhsEval>::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {};
+            try { val = FluidSystem::template fugacityCoefficient<FluidState, LhsEval>(fs, paramCache, phaseIdx, compIdx); } catch (...) {};
+            try { scalarVal = FluidSystem::template fugacityCoefficient<FluidState, Scalar>(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<decltype(tmpVal), RhsEval>::value, "The default return value must be the scalar used by the fluid state!"); } catch (...) {};
+            try { val = FluidSystem::template diffusionCoefficient<FluidState, LhsEval>(fs, paramCache, phaseIdx, compIdx); } catch (...) {};
+            try { scalarVal = FluidSystem::template fugacityCoefficient<FluidState, Scalar>(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);

tests/test_fluidsystems.cpp

@@ -24,6 +24,9 @@
  */
 #include "config.h"
 
+#include <opm/material/localad/Evaluation.hpp>
+#include <opm/material/localad/Math.hpp>
+
 #include "checkFluidSystem.hpp"
 
 // include all fluid systems in opm-material
@@ -45,8 +48,7 @@
 #include <opm/material/fluidstates/NonEquilibriumFluidState.hpp>
 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
 
-// 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 <opm/material/common/UniformTabulated2DFunction.hpp>
 
 namespace Opm {
@@ -77,116 +79,137 @@ public:
     };
 };
 
-int main(int argc, char **argv)
+// check the API of all fluid states
+template <class Scalar>
+void testAllFluidStates()
+{
+    typedef Opm::FluidSystems::H2ON2<Scalar, /*enableComplexRelations=*/false> FluidSystem;
+
+    // CompositionalFluidState
+    {   Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
+        checkFluidState<Scalar>(fs); }
+
+    // NonEquilibriumFluidState
+    {   Opm::NonEquilibriumFluidState<Scalar, FluidSystem> fs;
+        checkFluidState<Scalar>(fs); }
+
+    // ImmiscibleFluidState
+    {   Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
+        checkFluidState<Scalar>(fs); }
+
+    typedef Opm::CompositionalFluidState<Scalar, FluidSystem> BaseFluidState;
+    BaseFluidState baseFs;
+
+    // TemperatureOverlayFluidState
+    {   Opm::TemperatureOverlayFluidState<BaseFluidState> fs(baseFs);
+        checkFluidState<Scalar>(fs); }
+
+    // PressureOverlayFluidState
+    {   Opm::PressureOverlayFluidState<BaseFluidState> fs(baseFs);
+        checkFluidState<Scalar>(fs); }
+
+    // SaturationOverlayFluidState
+    {   Opm::SaturationOverlayFluidState<BaseFluidState> fs(baseFs);
+        checkFluidState<Scalar>(fs); }
+}
+
+template <class Scalar, class Evaluation, class LhsEval = Evaluation>
+void testAllFluidSystems()
 {
-    typedef double Scalar;
     typedef Opm::H2O<Scalar> H2O;
     typedef Opm::N2<Scalar> N2;
 
     typedef Opm::LiquidPhase<Scalar, H2O> Liquid;
     typedef Opm::GasPhase<Scalar, N2> Gas;
 
-    MyMpiHelper mpiHelper(argc, argv);
-
-    // check all fluid states
-    {
-        typedef Opm::FluidSystems::H2ON2<Scalar, /*enableComplexRelations=*/false> FluidSystem;
-
-        // CompositionalFluidState
-        {   Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
-            checkFluidState<Scalar>(fs); }
-
-        // NonEquilibriumFluidState
-        {   Opm::NonEquilibriumFluidState<Scalar, FluidSystem> fs;
-            checkFluidState<Scalar>(fs); }
-
-        // ImmiscibleFluidState
-        {   Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
-            checkFluidState<Scalar>(fs); }
-
-        typedef Opm::CompositionalFluidState<Scalar, FluidSystem> BaseFluidState;
-        BaseFluidState baseFs;
-
-        // TemperatureOverlayFluidState
-        {   Opm::TemperatureOverlayFluidState<BaseFluidState> fs(baseFs);
-            checkFluidState<Scalar>(fs); }
-
-        // PressureOverlayFluidState
-        {   Opm::PressureOverlayFluidState<BaseFluidState> fs(baseFs);
-            checkFluidState<Scalar>(fs); }
-
-        // SaturationOverlayFluidState
-        {   Opm::SaturationOverlayFluidState<BaseFluidState> fs(baseFs);
-            checkFluidState<Scalar>(fs); }
-    }
-
     // black-oil
-    {   typedef Opm::FluidSystems::BlackOil<Scalar> FluidSystem;
-        if (false) checkFluidSystem<Scalar, FluidSystem>(); }
+    {   typedef Opm::FluidSystems::BlackOil<Scalar, Evaluation> FluidSystem;
+        if (false) checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // Brine -- CO2
     {   typedef Opm::FluidSystems::BrineCO2<Scalar, Opm::FluidSystemsTest::CO2Tables> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // H2O -- N2
     {   typedef Opm::FluidSystems::H2ON2<Scalar, /*enableComplexRelations=*/false> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {   typedef Opm::FluidSystems::H2ON2<Scalar, /*enableComplexRelations=*/true> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // H2O -- N2 -- liquid phase
     {   typedef Opm::FluidSystems::H2ON2LiquidPhase<Scalar, /*enableComplexRelations=*/false> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {   typedef Opm::FluidSystems::H2ON2LiquidPhase<Scalar, /*enableComplexRelations=*/true> FluidSystem;
-         checkFluidSystem<Scalar, FluidSystem>(); }
+         checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // H2O -- Air
     {   typedef Opm::SimpleH2O<Scalar> H2O;
         const bool enableComplexRelations=false;
         typedef Opm::FluidSystems::H2OAir<Scalar, H2O, enableComplexRelations> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {   typedef Opm::SimpleH2O<Scalar> H2O;
         const bool enableComplexRelations=true;
         typedef Opm::FluidSystems::H2OAir<Scalar, H2O, enableComplexRelations> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {   typedef Opm::H2O<Scalar> H2O;
         const bool enableComplexRelations=false;
         typedef Opm::FluidSystems::H2OAir<Scalar, H2O, enableComplexRelations> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {   typedef Opm::H2O<Scalar> H2O;
         const bool enableComplexRelations=true;
         typedef Opm::FluidSystems::H2OAir<Scalar, H2O, enableComplexRelations> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // H2O -- Air -- Mesitylene
     {   typedef Opm::FluidSystems::H2OAirMesitylene<Scalar> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // H2O -- Air -- Xylene
     {   typedef Opm::FluidSystems::H2OAirXylene<Scalar> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // 2p-immiscible
     {   typedef Opm::FluidSystems::TwoPhaseImmiscible<Scalar, Liquid, Liquid> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {   typedef Opm::FluidSystems::TwoPhaseImmiscible<Scalar, Liquid, Gas> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {  typedef Opm::FluidSystems::TwoPhaseImmiscible<Scalar, Gas, Liquid> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     // 1p
     {   typedef Opm::FluidSystems::SinglePhase<Scalar, Liquid> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
 
     {   typedef Opm::FluidSystems::SinglePhase<Scalar, Gas> FluidSystem;
-        checkFluidSystem<Scalar, FluidSystem>(); }
+        checkFluidSystem<Scalar, FluidSystem, Evaluation, LhsEval>(); }
+}
+
+class TestAdTag;
+
+int main(int argc, char **argv)
+{
+    typedef double Scalar;
+    typedef Opm::LocalAd::Evaluation<Scalar, TestAdTag, 3> Evaluation;
+
+    MyMpiHelper mpiHelper(argc, argv);
+
+    // ensure that all fluid states are API-compliant
+    testAllFluidStates<Scalar>();
+    testAllFluidStates<Evaluation>();
+
+    // 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<Scalar, Scalar>();
+    testAllFluidSystems<Scalar, Evaluation>();
+    testAllFluidSystems<Scalar, Evaluation, Scalar>();
 
     return 0;
 }