Changes to add code for mixing in thermal simulation, for now setting…

… dones not change result

opm/material/common/UniformTabulated2DFunction.hpp

@@ -195,7 +195,7 @@ class UniformTabulated2DFunction
     Evaluation eval(const Evaluation& x, const Evaluation& y, bool extrapolate) const
     {
 #ifndef NDEBUG
-        if (!extrapolate && !applies(x,y)) {
+        if (!applies(x,y)) {
             std::string msg = "Attempt to get tabulated value for ("
                 +std::to_string(double(scalarValue(x)))+", "+std::to_string(double(scalarValue(y)))
                 +") on a table of extent "
@@ -212,6 +212,15 @@ class UniformTabulated2DFunction
             }
 
         };
+#else
+        if (!extrapolate && !applies(x,y)){
+            std::string msg = "Attempt to get tabulated value for ("
+                +std::to_string(double(scalarValue(x)))+", "+std::to_string(double(scalarValue(y)))
+                +") on a table of extent "
+                +std::to_string(xMin())+" to "+std::to_string(xMax())+" times "
+                +std::to_string(yMin())+" to "+std::to_string(yMax());
+            throw NumericalProblem(msg);
+        }
 #endif
 
         Evaluation alpha = xToI(x);

opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp

@@ -84,8 +84,8 @@ class BrineCo2Pvt
         : salinity_(salinity)
     {
         // Throw an error if reference state is not (T, p) = (15.56 C, 1 atm) = (288.71 K, 1.01325e5 Pa)
-        if (T_ref != Scalar(288.71) || P_ref != Scalar(1.01325e5)) {    
-            OPM_THROW(std::runtime_error, 
+        if (T_ref != Scalar(288.71) || P_ref != Scalar(1.01325e5)) {
+            OPM_THROW(std::runtime_error,
                 "BrineCo2Pvt class can only be used with default reference state (T, P) = (288.71 K, 1.01325e5 Pa)!");
         }
         setActivityModelSalt(activityModel);
@@ -174,6 +174,10 @@ class BrineCo2Pvt
     unsigned numRegions() const
     { return brineReferenceDensity_.size(); }
 
+    Scalar hVap(unsigned ) const{
+        return 0;
+    }
+
     /*!
      * \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
      */
@@ -314,8 +318,8 @@ class BrineCo2Pvt
                                                      const Evaluation& temperature,
                                                      const Evaluation& pressure) const
     {
- 	OPM_TIMEFUNCTION_LOCAL();
-	Evaluation rs_sat = rsSat(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
+    OPM_TIMEFUNCTION_LOCAL();
+    Evaluation rs_sat = rsSat(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
         return (1.0 - convertRsToXoG_(rs_sat,regionIdx)) * density(regionIdx, temperature, pressure, rs_sat, Evaluation(salinity_[regionIdx]))/brineReferenceDensity_[regionIdx];
     }
 
@@ -423,8 +427,8 @@ class BrineCo2Pvt
                        const Evaluation& Rs,
                        const Evaluation& salinity) const
     {
-	OPM_TIMEFUNCTION_LOCAL();        
-	Evaluation xlCO2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx), salinity);
+    OPM_TIMEFUNCTION_LOCAL();
+    Evaluation xlCO2 = convertXoGToxoG_(convertRsToXoG_(Rs,regionIdx), salinity);
         Evaluation result = liquidDensity_(temperature,
                                         pressure,
                                         xlCO2,
@@ -440,12 +444,12 @@ class BrineCo2Pvt
                      const Evaluation& pressure,
                      const Evaluation& salinity) const
     {
-	    OPM_TIMEFUNCTION_LOCAL();
+        OPM_TIMEFUNCTION_LOCAL();
         if (!enableDissolution_)
             return 0.0;
 
         // calulate the equilibrium composition for the given
-        // temperature and pressure. 
+        // temperature and pressure.
         Evaluation xgH2O;
         Evaluation xlCO2;
         BinaryCoeffBrineCO2::calculateMoleFractions(temperature,

opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.hpp

@@ -141,6 +141,10 @@ class BrineH2Pvt
     /*!
     * \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
     */
+    Scalar hVap(unsigned /*regionIdx*/) const{
+        return 0.0;
+    }
+
     template <class Evaluation>
     Evaluation internalEnergy(unsigned regionIdx,
                         const Evaluation& temperature,
@@ -236,7 +240,7 @@ class BrineH2Pvt
     {
         const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltconcentration);
         Evaluation rsSat = rsSat_(regionIdx, temperature, pressure, salinity);
-        return (1.0 - convertRsToXoG_(rsSat,regionIdx)) * 
+        return (1.0 - convertRsToXoG_(rsSat,regionIdx)) *
             density_(regionIdx, temperature, pressure, rsSat, salinity) / brineReferenceDensity_[regionIdx];
     }
 
@@ -251,7 +255,7 @@ class BrineH2Pvt
                                             const Evaluation& saltConcentration) const
     {
         const Evaluation salinity = salinityFromConcentration(regionIdx, temperature, pressure, saltConcentration);
-        return (1.0 - convertRsToXoG_(Rs,regionIdx)) * 
+        return (1.0 - convertRsToXoG_(Rs,regionIdx)) *
             density_(regionIdx, temperature, pressure, Rs, salinity) / brineReferenceDensity_[regionIdx];
     }
 
@@ -264,8 +268,8 @@ class BrineH2Pvt
                                             const Evaluation& pressure,
                                             const Evaluation& Rs) const
     {
-        return (1.0 - convertRsToXoG_(Rs, regionIdx)) * 
-            density_(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx])) / 
+        return (1.0 - convertRsToXoG_(Rs, regionIdx)) *
+            density_(regionIdx, temperature, pressure, Rs, Evaluation(salinity_[regionIdx])) /
                 brineReferenceDensity_[regionIdx];
     }
 
@@ -278,8 +282,8 @@ class BrineH2Pvt
                                                      const Evaluation& pressure) const
     {
         Evaluation rsSat = rsSat_(regionIdx, temperature, pressure, Evaluation(salinity_[regionIdx]));
-        return (1.0 - convertRsToXoG_(rsSat, regionIdx)) * 
-            density_(regionIdx, temperature, pressure, rsSat, Evaluation(salinity_[regionIdx])) / 
+        return (1.0 - convertRsToXoG_(rsSat, regionIdx)) *
+            density_(regionIdx, temperature, pressure, rsSat, Evaluation(salinity_[regionIdx])) /
                 brineReferenceDensity_[regionIdx];
     }
 
@@ -398,7 +402,7 @@ class BrineH2Pvt
 
     /*!
     * \brief Calculate density of aqueous solution (H2O-NaCl/brine and H2).
-    * 
+    *
     * \param temperature temperature [K]
     * \param pressure pressure [Pa]
     * \param Rs gas dissolution factor [-]
@@ -426,7 +430,7 @@ class BrineH2Pvt
     /*!
     * \brief Calculated the density of the aqueous solution where contributions of salinity and dissolved H2 is taken
     * into account.
-    * 
+    *
     * \param T temperature [K]
     * \param pl liquid pressure [Pa]
     * \param xlH2 mole fraction H2 [-]
@@ -444,9 +448,9 @@ class BrineH2Pvt
 
         // check if pressure and temperature is valid
         if(!extrapolate && T < 273.15) {
-            const std::string msg = 
+            const std::string msg =
                 "Liquid density for Brine and H2 is only "
-                "defined above 273.15K (is " + 
+                "defined above 273.15K (is " +
                 std::to_string(getValue(T)) + "K)";
             throw NumericalProblem(msg);
         }
@@ -470,7 +474,7 @@ class BrineH2Pvt
     /*!
     * \brief Density of aqueous solution with dissolved H2. Formula from Li et al. (2018) and Garica, Lawrence Berkeley
     * National Laboratory, 2001.
-    * 
+    *
     * \param temperature [K]
     * \param pl liquid pressure [Pa]
     * \param xlH2 mole fraction [-]
@@ -503,7 +507,7 @@ class BrineH2Pvt
     /*!
     * \brief Convert a gas dissolution factor to the the corresponding mass fraction of the gas component in the oil
     * phase.
-    * 
+    *
     * \param Rs gass dissolution factor [-]
     * \param regionIdx region index
     */
@@ -547,7 +551,7 @@ class BrineH2Pvt
     /*!
     * \brief Convert the mass fraction of the gas component in the oil phase to the corresponding gas dissolution
     * factor.
-    * 
+    *
     * \param XoG mass fraction [-]
     * \param regionIdx region index
     */
@@ -579,7 +583,7 @@ class BrineH2Pvt
 
         // calulate the equilibrium composition for the given temperature and pressure
         LhsEval xlH2 = BinaryCoeffBrineH2::calculateMoleFractions(temperature, pressure, salinity, extrapolate);
-        
+
         // normalize the phase compositions
         xlH2 = max(0.0, min(1.0, xlH2));
 
@@ -654,9 +658,9 @@ class BrineH2Pvt
     }
 
     template <class LhsEval>
-    const LhsEval salinityFromConcentration(unsigned regionIdx, 
-                                            const LhsEval&T, 
-                                            const LhsEval& P, 
+    const LhsEval salinityFromConcentration(unsigned regionIdx,
+                                            const LhsEval&T,
+                                            const LhsEval& P,
                                             const LhsEval& saltConcentration) const
     {
         if (enableSaltConcentration_)

opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp

@@ -71,7 +71,7 @@ class Co2GasPvt
     {
         // Throw an error if reference state is not (T, p) = (15.56 C, 1 atm) = (288.71 K, 1.01325e5 Pa)
         if (T_ref != Scalar(288.71) || P_ref != Scalar(1.01325e5)) {
-            OPM_THROW(std::runtime_error, 
+            OPM_THROW(std::runtime_error,
                 "BrineCo2Pvt class can only be used with default reference state (T, P) = (288.71 K, 1.01325e5 Pa)!");
         }
         setActivityModelSalt(activityModel);
@@ -146,6 +146,9 @@ class Co2GasPvt
     unsigned numRegions() const
     { return gasReferenceDensity_.size(); }
 
+    Scalar hVap(unsigned ) const{
+        return 0;
+    }
     /*!
      * \brief Returns the specific enthalpy [J/kg] of gas given a set of parameters.
      */
@@ -158,7 +161,7 @@ class Co2GasPvt
     {
         OPM_TIMEBLOCK_LOCAL(internalEnergy);
         // use the gasInternalEnergy of CO2
-        return CO2::gasInternalEnergy(temperature, pressure, extrapolate); 
+        return CO2::gasInternalEnergy(temperature, pressure, extrapolate);
 
         // account for H2O in the gas phase
         // Evaluation result = 0;

opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityBrinePvt.hpp

@@ -107,6 +107,9 @@ class ConstantCompressibilityBrinePvt
         throw std::runtime_error("Requested the enthalpy of water but the thermal option is not enabled");
     }
 
+    Scalar hVap(unsigned) const{
+        throw std::runtime_error("Requested the hvap of oil but the thermal option is not enabled");
+    }
     /*!
      * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
      */

opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityOilPvt.hpp

@@ -140,7 +140,9 @@ class ConstantCompressibilityOilPvt
     {
         throw std::runtime_error("Requested the enthalpy of oil but the thermal option is not enabled");
     }
-
+    Scalar hVap(unsigned) const{
+        throw std::runtime_error("Requested the hvap of oil but the thermal option is not enabled");
+    }
     /*!
      * \brief Returns the dynamic viscosity [Pa s] of gas saturated oil given a pressure
      *        and a phase composition.

opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp

@@ -140,7 +140,9 @@ class ConstantCompressibilityWaterPvt
     {
         throw std::runtime_error("Requested the enthalpy of water but the thermal option is not enabled");
     }
-
+    Scalar hVap(unsigned) const{
+        throw std::runtime_error("Requested the hvap of oil but the thermal option is not enabled");
+    }
 
     /*!
      * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
@@ -232,9 +234,9 @@ class ConstantCompressibilityWaterPvt
 
         // TODO (?): consider the salt concentration of the brine
         bw = (1.0 + X*(1.0 + X/2.0))/BwRef;
-        
+
         Scalar BwMuwRef = waterViscosity_[regionIdx]*BwRef;
-        
+
         const Evaluation& Y =
             (waterCompressibility_[regionIdx] - waterViscosibility_[regionIdx])
             * (pressure - pRef);
@@ -261,7 +263,7 @@ class ConstantCompressibilityWaterPvt
     {
         throw std::runtime_error("Not implemented: The PVT model does not provide a diffusionCoefficient()");
     }
-    
+
     /*!
      * \brief Returns the gas dissolution factor \f$R_s\f$ [m^3/m^3] of the water phase.
      */

opm/material/fluidsystems/blackoilpvt/DeadOilPvt.hpp

@@ -114,6 +114,9 @@ class DeadOilPvt
         throw std::runtime_error("Requested the enthalpy of oil but the thermal option is not enabled");
     }
 
+    Scalar hVap(unsigned) const{
+        throw std::runtime_error("Requested the hvap of oil but the thermal option is not enabled");
+    }
     /*!
      * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
      */

opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp

@@ -126,7 +126,9 @@ class DryGasPvt
     {
         throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
     }
-
+    Scalar hVap(unsigned) const{
+        throw std::runtime_error("Requested the hvap of oil but the thermal option is not enabled");
+    }
     /*!
      * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
      */
@@ -199,10 +201,10 @@ class DryGasPvt
     template <class Evaluation = Scalar>
     Evaluation saturatedWaterVaporizationFactor(unsigned /*regionIdx*/,
                                               const Evaluation& /*temperature*/,
-                                              const Evaluation& /*pressure*/, 
+                                              const Evaluation& /*pressure*/,
                                               const Evaluation& /*saltConcentration*/) const
     { return 0.0; }
-    
+
 
     /*!
      * \brief Returns the oil vaporization factor \f$R_v\f$ [m^3/m^3] of the oil phase.

opm/material/fluidsystems/blackoilpvt/DryHumidGasPvt.hpp

@@ -151,7 +151,9 @@ class DryHumidGasPvt
     {
         throw std::runtime_error("Requested the enthalpy of gas but the thermal option is not enabled");
     }
-
+    Scalar hVap(unsigned) const{
+        throw std::runtime_error("Requested the hvap of oil but the thermal option is not enabled");
+    }
     /*!
      * \brief Returns the dynamic viscosity [Pa s] of the fluid phase given a set of parameters.
      */