Speedup the temperature module

- Use its own fluidstate and Evaluation class to avoid having to add an extra derivative to the main equations
- Sotre the scaling factor to avoid calling the parameter all the time

flow/flow_biofilm.cpp

@@ -91,13 +91,11 @@ private:
     // get the energy module type to determine the equation number
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_blackoil_temp.cpp

@@ -49,8 +49,7 @@ template<class TypeTag>
 struct EnableDiffusion<TypeTag, TTag::flowBlackoilTempProblem> { static constexpr bool value = false; };
 
 template<class TypeTag>
-struct NumDerivatives<TypeTag, TTag::flowBlackoilTempProblem>
-{ static constexpr int value = GetPropType<TypeTag, Properties::Indices>::numDerivatives; };
+struct AvoidElementContext<TypeTag, TTag::flowBlackoilTempProblem> { static constexpr bool value = true; };
 
 
 }}

flow/flow_gasoil.cpp

@@ -62,13 +62,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
   using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                        getPropValue<TypeTag, Properties::EnableExtbo>(),
                                        getPropValue<TypeTag, Properties::EnablePolymer>(),
                                        numEnergyVars,
-                                       enableSeqImpEnergy,
                                        getPropValue<TypeTag, Properties::EnableFoam>(),
                                        getPropValue<TypeTag, Properties::EnableBrine>(),
                                        /*PVOffset=*/0,

flow/flow_gasoil_energy.cpp

@@ -49,14 +49,12 @@ private:
     // get the energy module type to determine the equation number
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 
 public:
   using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                        getPropValue<TypeTag, Properties::EnableExtbo>(),
                                        getPropValue<TypeTag, Properties::EnablePolymer>(),
                                        numEnergyVars,
-                                       enableSeqImpEnergy,
                                        getPropValue<TypeTag, Properties::EnableFoam>(),
                                        getPropValue<TypeTag, Properties::EnableBrine>(),
                                        /*PVOffset=*/0,

flow/flow_gasoildiffuse.cpp

@@ -64,13 +64,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
   using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                        getPropValue<TypeTag, Properties::EnableExtbo>(),
                                        getPropValue<TypeTag, Properties::EnablePolymer>(),
                                        numEnergyVars,
-                                       enableSeqImpEnergy,
                                        getPropValue<TypeTag, Properties::EnableFoam>(),
                                        getPropValue<TypeTag, Properties::EnableBrine>(),
                                        /*PVOffset=*/0,

flow/flow_gaswater.cpp

@@ -59,13 +59,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_gaswater_brine.cpp

@@ -77,13 +77,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_gaswater_dissolution.cpp

@@ -75,13 +75,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_gaswater_dissolution_diffuse.cpp

@@ -78,13 +78,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_gaswater_energy.cpp

@@ -78,13 +78,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_gaswater_saltprec_energy.cpp

@@ -80,13 +80,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_gaswater_saltprec_vapwat.cpp

@@ -80,13 +80,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_gaswater_solvent.cpp

@@ -57,13 +57,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_micp.cpp

@@ -60,13 +60,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_oilwater.cpp

@@ -58,14 +58,12 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_oilwater_brine.cpp

@@ -62,13 +62,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_oilwater_polymer.cpp

@@ -48,13 +48,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_oilwater_polymer_injectivity.cpp

@@ -55,13 +55,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilTwoPhaseIndices<0,
                                          0,
                                          2,
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_onephase.cpp

@@ -55,13 +55,11 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 public:
     using type = BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                          getPropValue<TypeTag, Properties::EnableExtbo>(),
                                          getPropValue<TypeTag, Properties::EnablePolymer>(),
                                          numEnergyVars,
-                                         enableSeqImpEnergy,
                                          getPropValue<TypeTag, Properties::EnableFoam>(),
                                          getPropValue<TypeTag, Properties::EnableBrine>(),
                                          /*PVOffset=*/0,

flow/flow_onephase_energy.cpp

@@ -50,14 +50,12 @@ private:
     using FluidSystem = GetPropType<BaseTypeTag, Properties::FluidSystem>;
     static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
     static constexpr int numEnergyVars = energyModuleType == EnergyModules::FullyImplicitThermal;
-    static constexpr bool enableSeqImpEnergy = energyModuleType == EnergyModules::SequentialImplicitThermal;
 
 public:
     using type = Opm::BlackOilOnePhaseIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
                                               getPropValue<TypeTag, Properties::EnableExtbo>(),
                                               getPropValue<TypeTag, Properties::EnablePolymer>(),
                                               numEnergyVars,
-                                              enableSeqImpEnergy,
                                               getPropValue<TypeTag, Properties::EnableFoam>(),
                                               getPropValue<TypeTag, Properties::EnableBrine>(),
                                               /*PVOffset=*/0,

flowexperimental/BlackOilEnergyIntensiveQuantitiesGlobalIndex.hpp

@@ -117,17 +117,17 @@ class BlackOilEnergyIntensiveQuantitiesGlobalIndex<TypeTag, EnergyModules::Seque
     using ParamCache = typename FluidSystem::template ParameterCache<Evaluation>;
 
     using Indices = GetPropType<TypeTag, Properties::Indices>;
-    static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
     static constexpr unsigned numPhases = FluidSystem::numPhases;
 public:
+
     void updateTemperature_([[maybe_unused]] const Problem& problem,
                             [[maybe_unused]] const PrimaryVariables& priVars,
                             [[maybe_unused]] unsigned globalSpaceIndex,
                             [[maybe_unused]] unsigned timeIdx)
     {
-        auto& fs = this->asImp_().fluidState_;
-        Evaluation T = Evaluation::createVariable(problem.temperature(globalSpaceIndex, timeIdx), Indices::temperatureIdx);
-        fs.setTemperature(T);
+        throw std::logic_error("updateTemperature not implemented "
+                        "SequentialImplicitThermal can not be used with"
+                        "global intensive quantites yet.");
     }
 
     void updateEnergyQuantities_([[maybe_unused]] const Problem& problem,
@@ -136,32 +136,11 @@ public:
                                  [[maybe_unused]] unsigned timeIdx,
                                  const ParamCache& paramCache)
     {
-        auto& fs = Parent::asImp_().fluidState_;
-
-        // compute the specific enthalpy of the fluids, the specific enthalpy of the rock
-        // and the thermal conductivity coefficients
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-            if (!FluidSystem::phaseIsActive(phaseIdx)) {
-                continue;
-            }
-
-            const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
-            fs.setEnthalpy(phaseIdx, h);
-        }
-
-        const auto& solidEnergyLawParams = problem().solidEnergyLawParams(globalSpaceIndex, timeIdx);
-        this->rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
-
-        const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIndex, timeIdx);
-        this->totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
-
-        // Retrieve the rock fraction from the problem
-        // Usually 1 - porosity, but if pvmult is used to modify porosity
-        // we will apply the same multiplier to the rock fraction
-        // i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
-        // multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
-        this->rockFraction_ = problem.rockFraction(globalSpaceIndex, timeIdx);
+        throw std::logic_error("updateEnergyQuantities_ not implemented "
+                "SequentialImplicitThermal can not be used with"
+                "global intensive quantites yet.");
     }
+
 };
 
 template <class TypeTag>

flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp

@@ -133,8 +133,8 @@ class BlackOilIntensiveQuantitiesGlobalIndex
 public:
     using FluidState = BlackOilFluidState<Evaluation,
                                           FluidSystem,
-                                          energyModuleType == EnergyModules::ConstantTemperature,
-                                          (energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
+                                          energyModuleType != EnergyModules::NoTemperature,
+                                          energyModuleType == EnergyModules::FullyImplicitThermal,
                                           compositionSwitchEnabled,
                                           enableVapwat,
                                           enableBrine,

opm/models/blackoil/blackoilenergymodules.hh

@@ -32,6 +32,7 @@
 
 #include <opm/material/common/Tabulated1DFunction.hpp>
 #include <opm/material/common/Valgrind.hpp>
+#include <opm/material/fluidstates/BlackOilFluidState.hpp>
 
 #include <opm/models/blackoil/blackoilproperties.hh>
 #include <opm/models/common/quantitycallbacks.hh>
@@ -523,15 +524,14 @@ class BlackOilEnergyIntensiveQuantities<TypeTag, EnergyModules::SequentialImplic
     using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
-    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
 
 public:
 
     void updateTemperature_(const Problem& problem, unsigned globalDofIdx, unsigned timeIdx)
     {
+        // update the temperature for output (without derivatives)
         auto& fs = asImp_().fluidState_;
-        const Evaluation T = Evaluation::createVariable(problem.temperature(globalDofIdx, timeIdx), Indices::temperatureIdx);
-        fs.setTemperature(T);
+        fs.setTemperature(problem.temperature(globalDofIdx, timeIdx));
     }
 
     void updateTemperature_(const ElementContext& elemCtx,
@@ -559,56 +559,30 @@ public:
                                  unsigned dofIdx,
                                  unsigned timeIdx)
     {
-        updateEnergyQuantities_(elemCtx.problem(), elemCtx.globalSpaceIndex(dofIdx, timeIdx), timeIdx);
     }
 
     void updateEnergyQuantities_(const Problem& problem,
                                  const unsigned globalSpaceIdx,
                                  const unsigned timeIdx)
     {
-        auto& fs = asImp_().fluidState_;
-
-        // compute the specific enthalpy of the fluids, the specific enthalpy of the rock
-        // and the thermal conductivity coefficients
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
-            if (!FluidSystem::phaseIsActive(phaseIdx)) {
-                continue;
-            }
-
-            const auto& h = FluidSystem::enthalpy(fs, phaseIdx, problem.pvtRegionIndex(globalSpaceIdx));
-            fs.setEnthalpy(phaseIdx, h);
-        }
-
-        const auto& solidEnergyLawParams = problem.solidEnergyLawParams(globalSpaceIdx, timeIdx);
-        rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
-
-        const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIdx, timeIdx);
-        totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
-
-        // Retrieve the rock fraction from the problem
-        // Usually 1 - porosity, but if pvmult is used to modify porosity
-        // we will apply the same multiplier to the rock fraction
-        // i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
-        // multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
-        rockFraction_ = problem.rockFraction(globalSpaceIdx, timeIdx);
     }
 
     const Evaluation& rockInternalEnergy() const
-    { return rockInternalEnergy_; }
+    {
+        throw std::logic_error("Requested the rock internal energy, which is "
+                             "unavailable because energy is not conserved");
+    }
 
     const Evaluation& totalThermalConductivity() const
-    { return totalThermalConductivity_; }
-
-    const Scalar& rockFraction() const
-    { return rockFraction_; }
+    {
+        throw std::logic_error("Requested the total thermal conductivity, which is "
+                             "unavailable because energy is not conserved");
+    }
 
 protected:
     Implementation& asImp_()
     { return *static_cast<Implementation*>(this); }
 
-    Evaluation rockInternalEnergy_;
-    Evaluation totalThermalConductivity_;
-    Scalar rockFraction_;
 };
 
 template <class TypeTag>
@@ -678,10 +652,9 @@ class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::FullyImplicitThe
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
-    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
 
 public:
-    template<class FluidState>
+    template<class Evaluation, class FluidState, class IntensiveQuantities>
     static void updateEnergy(Evaluation& energyFlux,
                              const unsigned& focusDofIndex,
                              const unsigned& inIdx,
@@ -786,7 +759,7 @@ public:
         updateEnergyBoundary(energyFlux_, inIq, focusDofIdx, inIdx, alpha, boundaryFs);
     }
 
-    template <class BoundaryFluidState>
+    template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
     static void updateEnergyBoundary(Evaluation& energyFlux,
                                      const IntensiveQuantities& inIq,
                                      unsigned focusDofIndex,
@@ -839,12 +812,11 @@ template <class TypeTag>
 class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::ConstantTemperature>
 {
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
-    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
-    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
 
 public:
-    template<class FluidState>
+    template<class Evaluation, class FluidState, class IntensiveQuantities>
     static void updateEnergy(Evaluation& /*energyFlux*/,
                              const unsigned& /*focusDofIndex*/,
                              const unsigned& /*inIdx*/,
@@ -870,7 +842,7 @@ public:
                               const BoundaryFluidState&)
     {}
 
-    template <class BoundaryFluidState>
+    template <class BoundaryFluidState,class IntensiveQuantities>
     static void updateEnergyBoundary(Evaluation& /*heatFlux*/,
                                      const IntensiveQuantities& /*inIq*/,
                                      unsigned /*focusDofIndex*/,
@@ -888,12 +860,11 @@ template <class TypeTag>
 class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::SequentialImplicitThermal>
 {
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
-    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
-    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
 
 public:
-    template<class FluidState>
+    template<class Evaluation, class FluidState, class IntensiveQuantities>
     static void updateEnergy(Evaluation& energyFlux,
                              const unsigned& focusDofIndex,
                              const unsigned& inIdx,
@@ -932,7 +903,7 @@ public:
                               const BoundaryFluidState&)
     { }
 
-    template <class BoundaryFluidState>
+    template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
     static void updateEnergyBoundary(Evaluation& /*heatFlux*/,
                                      const IntensiveQuantities& /*inIq*/,
                                      unsigned /*focusDofIndex*/,
@@ -951,11 +922,10 @@ class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::NoTemperature>
 {
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
-    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
 public:
-    template<class FluidState>
+    template<class Evaluation, class FluidState, class IntensiveQuantities>
     static void updateEnergy(Evaluation& /*energyFlux*/,
                              const unsigned& /*focusDofIndex*/,
                              const unsigned& /*inIdx*/,
@@ -981,7 +951,7 @@ public:
                               const BoundaryFluidState&)
     {}
 
-    template <class BoundaryFluidState>
+    template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
     static void updateEnergyBoundary(Evaluation& /*heatFlux*/,
                                      const IntensiveQuantities& /*inIq*/,
                                      unsigned /*focusDofIndex*/,

opm/models/blackoil/blackoilintensivequantities.hh

@@ -137,8 +137,8 @@ class BlackOilIntensiveQuantities
 public:
     using FluidState = BlackOilFluidState<Evaluation,
                                           FluidSystem,
-                                          energyModuleType == EnergyModules::ConstantTemperature,
-                                          (energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
+                                          energyModuleType != EnergyModules::NoTemperature,
+                                          energyModuleType == EnergyModules::FullyImplicitThermal,
                                           compositionSwitchEnabled,
                                           enableVapwat,
                                           enableBrine,
@@ -147,9 +147,8 @@ public:
                                           Indices::numPhases>;
     using ScalarFluidState = BlackOilFluidState<Scalar,
                                                 FluidSystem,
-                                                energyModuleType == EnergyModules::ConstantTemperature,
-                                                (energyModuleType == EnergyModules::FullyImplicitThermal || 
-                                                    energyModuleType == EnergyModules::SequentialImplicitThermal),
+                                                energyModuleType != EnergyModules::NoTemperature,
+                                                energyModuleType == EnergyModules::FullyImplicitThermal,
                                                 compositionSwitchEnabled,
                                                 enableVapwat,
                                                 enableBrine,
@@ -633,8 +632,7 @@ public:
         if constexpr (enablePolymer) {
             asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
         }
-        if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal ||
-                      energyModuleType == EnergyModules::SequentialImplicitThermal) {
+        if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
             asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx);
         }
         if constexpr (enableFoam) {

opm/models/blackoil/blackoilmodel.hh

@@ -145,7 +145,6 @@ struct Indices<TypeTag, TTag::BlackOilModel>
                                                     getPropValue<TypeTag, Properties::EnableExtbo>(),
                                                     getPropValue<TypeTag, Properties::EnablePolymer>(),
                                                     getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal,
-                                                    getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal,
                                                     getPropValue<TypeTag, Properties::EnableFoam>(),
                                                     getPropValue<TypeTag, Properties::EnableBrine>(),
                                                     /*PVOffset=*/0,

opm/models/blackoil/blackoilonephaseindices.hh

@@ -43,7 +43,6 @@ template<unsigned numSolventsV,
          unsigned numExtbosV,
          unsigned numPolymersV,
          unsigned numEnergyV,
-         bool enableSequentialImplicitThermal,
          bool enableFoam,
          bool enableBrine,
          unsigned PVOffset,
@@ -105,8 +104,8 @@ struct BlackOilOnePhaseIndices
     static constexpr int numEq = numPhases + numSolvents + numExtbos + numPolymers +
                                  numEnergy + numFoam + numBrine + numMICPs;
 
-    //! For seqential implicit approches we evaluate the intensive quantities with a larger number of derivatives than equations
-    static constexpr int numDerivatives = numEq + enableSequentialImplicitThermal;
+    //! We may want to evaluate with more derivatives that equations, but not this time
+    static constexpr int numDerivatives = numEq;
 
     //////////////////////////////
     // Primary variable indices
@@ -184,7 +183,7 @@ struct BlackOilOnePhaseIndices
 
     //! Index of the primary variable for temperature
     static constexpr int temperatureIdx  =
-        (enableSequentialImplicitThermal || enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000;
+        (enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000;
 
     //////////////////////
     // Equation indices

opm/models/blackoil/blackoiltwophaseindices.hh

@@ -43,7 +43,6 @@ template<unsigned numSolventsV,
          unsigned numExtbosV,
          unsigned numPolymersV,
          unsigned numEnergyV,
-         bool enableSequentialImplicitThermal,
          bool enableFoam,
          bool enableBrine,
          int PVOffset,
@@ -105,8 +104,8 @@ struct BlackOilTwoPhaseIndices
     static constexpr int numEq = numPhases + numSolvents + numExtbos + numPolymers +
                                  numEnergy + numFoam + numBrine + numBioComp;
 
-    //! For seqential implicit approches we evaluate the intensive quantities with a larger number of derivatives than equations
-    static constexpr int numDerivatives = numEq + enableSequentialImplicitThermal;
+    //! We may want to evaluate with more derivatives that equations, but not this time
+    static constexpr int numDerivatives = numEq;
 
     //////////////////////////////
     // Primary variable indices
@@ -177,7 +176,7 @@ struct BlackOilTwoPhaseIndices
 
     //! Index of the primary variable for temperature
     static constexpr int temperatureIdx  =
-        (enableSequentialImplicitThermal || enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numBioComp + numFoam + numBrine : - 1000;
+        (enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numBioComp + numFoam + numBrine : - 1000;
 
     //////////////////////
     // Equation indices

opm/models/blackoil/blackoilvariableandequationindices.hh

@@ -39,7 +39,6 @@ template<unsigned numSolventsV,
          unsigned numExtbosV,
          unsigned numPolymersV,
          unsigned numEnergyV,
-         bool enableSequentialImplicitThermal,
          bool enableFoam,
          bool enableBrine,
          int PVOffset,
@@ -93,8 +92,8 @@ struct BlackOilVariableAndEquationIndices
     static constexpr int numEq = numPhases + numSolvents + numExtbos + numPolymers +
                                  numEnergy + numFoam + numBrine;
 
-    //! For seqential implicit approches we evaluate the intensive quantities with a larger number of derivatives than equations
-    static constexpr int numDerivatives = numEq + enableSequentialImplicitThermal;
+    //! We may want to evaluate with more derivatives that equations, but not this time
+    static constexpr int numDerivatives = numEq;
 
 
     ////////
@@ -160,7 +159,7 @@ struct BlackOilVariableAndEquationIndices
 
     //! Index of the primary variable for temperature
     static constexpr int temperatureIdx  =
-        (enableSequentialImplicitThermal || enableFullyImplicitThermal) ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
+        enableFullyImplicitThermal ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
 
 
     ////////

opm/simulators/aquifers/AquiferAnalytical.hpp

@@ -78,8 +78,8 @@ public:
 
     using FluidState = BlackOilFluidState<Eval,
                                           FluidSystem,
-                                          energyModuleType == EnergyModules::ConstantTemperature,
-                                          (energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
+                                          energyModuleType != EnergyModules::NoTemperature,
+                                          energyModuleType == EnergyModules::FullyImplicitThermal,
                                           BlackoilIndices::gasEnabled,
                                           enableVapwat,
                                           enableBrine,

opm/simulators/flow/EquilInitializer.hpp

@@ -88,8 +88,8 @@ public:
     // internal energy!
     using ScalarFluidState = BlackOilFluidState<Scalar,
                                                 FluidSystem,
-                                                energyModuleType == EnergyModules::ConstantTemperature,
-                                                (energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
+                                                energyModuleType != EnergyModules::NoTemperature,
+                                                energyModuleType == EnergyModules::FullyImplicitThermal,
                                                 enableDissolution,
                                                 enableVapwat,
                                                 enableBrine,
@@ -173,7 +173,7 @@ public:
                 const auto& rho = FluidSystem::density(fluidState, phaseIdx, regionIdx);
                 fluidState.setDensity(phaseIdx, rho);
 
-                if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal) {
+                if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
                     const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, regionIdx);
                     fluidState.setEnthalpy(phaseIdx, h);
                 }

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -857,7 +857,7 @@ public:
 
                     const auto& rho = FluidSystem::density(fluidState, phaseIdx, pvtRegionIdx);
                     fluidState.setDensity(phaseIdx, rho);
-                    if constexpr (energyModuleType == EnergyModules::SequentialImplicitThermal || energyModuleType == EnergyModules::FullyImplicitThermal) {
+                    if constexpr (energyModuleType == EnergyModules::FullyImplicitThermal) {
                         const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, pvtRegionIdx);
                         fluidState.setEnthalpy(phaseIdx, h);
                     }

opm/simulators/flow/GenericTemperatureModel_impl.hpp

@@ -115,8 +115,11 @@ void GenericTemperatureModel<Grid,GridView,DofMapper,Stencil,FluidSystem,Scalar>
 doInit(std::size_t numGridDof)
 {
     doTemp_ = eclState_.getSimulationConfig().isTemp();
-
     temperature_.resize(numGridDof);
+
+    if (!doTemp_)
+        return;
+
     energyVector_.resize(numGridDof);
     // allocate matrix for storing the Jacobian of the temperature residual
     energyMatrix_ = std::make_unique<EnergyMatrix>(numGridDof, numGridDof, EnergyMatrix::random);

opm/simulators/flow/TemperatureModel.hpp

@@ -38,6 +38,7 @@
 #include <opm/simulators/linalg/findOverlapRowsAndColumns.hpp>
 #include <opm/simulators/aquifers/AquiferGridUtils.hpp>
 
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -58,6 +59,108 @@ namespace Opm {
 
 template<typename Scalar, typename IndexTraits> class WellState;
 
+
+template <class TypeTag>
+class BlackOilEnergyIntensiveQuantitiesTemp
+{
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
+    using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+
+    using EvaluationTemp = DenseAd::Evaluation<Scalar,1>;
+
+
+    enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
+    enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
+    enum { enableDisgasInWater = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
+    enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
+    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
+    static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
+    using FluidStateTemp = BlackOilFluidState<  EvaluationTemp,
+                                                FluidSystem,
+                                                true, //store temperature
+                                                true, // store enthalpy
+                                                compositionSwitchEnabled,
+                                                enableVapwat,
+                                                enableBrine,
+                                                enableSaltPrecipitation,
+                                                enableDisgasInWater,
+                                                Indices::numPhases>;
+
+public:
+
+    void updateTemperature_(const Problem& problem, unsigned globalDofIdx, unsigned timeIdx)
+    {
+        const EvaluationTemp T = EvaluationTemp::createVariable(problem.temperature(globalDofIdx, timeIdx), 0);
+        fluidState_.setTemperature(T);
+    }
+
+    void updateEnergyQuantities_(const Problem& problem,
+                                 const unsigned globalSpaceIdx,
+                                 const unsigned timeIdx)
+    {
+        // compute the specific enthalpy of the fluids, the specific enthalpy of the rock
+        // and the thermal conductivity coefficients
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx)) {
+                continue;
+            }
+
+            const auto& h = FluidSystem::enthalpy(fluidState_, phaseIdx, problem.pvtRegionIndex(globalSpaceIdx));
+            fluidState_.setEnthalpy(phaseIdx, h);
+        }
+
+        const auto& solidEnergyLawParams = problem.solidEnergyLawParams(globalSpaceIdx, timeIdx);
+        rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fluidState_);
+
+        const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIdx, timeIdx);
+        totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fluidState_);
+
+        // Retrieve the rock fraction from the problem
+        // Usually 1 - porosity, but if pvmult is used to modify porosity
+        // we will apply the same multiplier to the rock fraction
+        // i.e. pvmult*(1 - porosity) and thus interpret multpv as a volume
+        // multiplier. This is to avoid negative rock volume for pvmult*porosity > 1
+        rockFraction_ = problem.rockFraction(globalSpaceIdx, timeIdx);
+    }
+
+    const EvaluationTemp& rockInternalEnergy() const
+    { return rockInternalEnergy_; }
+
+    const EvaluationTemp& totalThermalConductivity() const
+    { return totalThermalConductivity_; }
+
+    const Scalar& rockFraction() const
+    { return rockFraction_; }
+
+    const FluidStateTemp& fluidStateTemp() const
+    { return fluidState_; }
+
+    template <class FluidState>
+    void setFluidState(const FluidState& fs)
+    {
+        // copy the needed part of the fluid state
+        fluidState_.setPvtRegionIndex(fs.pvtRegionIndex());
+        fluidState_.setRs(getValue(fs.Rs()));
+        fluidState_.setRv(getValue(fs.Rv()));
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+            fluidState_.setPressure(phaseIdx, getValue(fs.pressure(phaseIdx)));
+            fluidState_.setDensity(phaseIdx, getValue(fs.density(phaseIdx)));
+            fluidState_.setSaturation(phaseIdx, getValue(fs.saturation(phaseIdx)));
+            fluidState_.setInvB(phaseIdx, getValue(fs.invB(phaseIdx)));
+        }
+    }
+
+protected:
+    EvaluationTemp rockInternalEnergy_;
+    EvaluationTemp totalThermalConductivity_;
+    FluidStateTemp fluidState_;
+    Scalar rockFraction_;
+};
+
 /*!
  * \ingroup BlackOilSimulator
  *
@@ -86,14 +189,15 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
-    using TemperatureEvaluation = DenseAd::Evaluation<Scalar,1>;
-    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using Evaluation = DenseAd::Evaluation<Scalar,1>;
     using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
     using IndexTraits = typename FluidSystem::IndexTraitsType;
     using WellStateType = WellState<Scalar, IndexTraits>;
 
+    using IntensiveQuantitiesTemp = BlackOilEnergyIntensiveQuantitiesTemp<TypeTag>;
+
     using EnergyMatrix = typename BaseType::EnergyMatrix;
     using EnergyVector = typename BaseType::EnergyVector;
 
@@ -102,6 +206,7 @@ class TemperatureModel : public GenericTemperatureModel<GetPropType<TypeTag, Pro
     enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
     enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    static constexpr int temperatureIdx = 0;
 
 public:
     explicit TemperatureModel(Simulator& simulator)
@@ -134,6 +239,9 @@ public:
         // find and store the overlap cells
         const auto& elemMapper = simulator_.model().elementMapper();
         detail::findOverlapAndInterior(simulator_.vanguard().grid(), elemMapper, overlapRows_, interiorRows_);
+
+        // set the scaling factor
+        scalingFactor_ = getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
     void beginTimeStep()
@@ -142,10 +250,13 @@ public:
             return;
         }
 
-        // We copy the intensive quantities here to make it possible to update them
+        // We use the specialized intensive quantities here with only the temperature derivative
         const unsigned int numCells = simulator_.model().numTotalDof();
+        #ifdef _OPENMP
+        #pragma omp parallel for
+        #endif
         for (unsigned globI = 0; globI < numCells; ++globI) {
-            intQuants_[globI] = simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0);
+            intQuants_[globI].setFluidState(simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0).fluidState());
             intQuants_[globI].updateTemperature_(simulator_.problem(), globI, /*timeIdx*/ 0);
             intQuants_[globI].updateEnergyQuantities_(simulator_.problem(), globI, /*timeIdx*/ 0);
         }
@@ -164,10 +275,13 @@ public:
             return;
         }
 
-        // We copy the intensive quantities here to make it possible to update them
+        // We use the specialized intensive quantities here with only the temperature derivative
         const unsigned int numCells = simulator_.model().numTotalDof();
+        #ifdef _OPENMP
+        #pragma omp parallel for
+        #endif
         for (unsigned globI = 0; globI < numCells; ++globI) {
-            intQuants_[globI] = simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0);
+            intQuants_[globI].setFluidState(simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0).fluidState());
             intQuants_[globI].updateTemperature_(simulator_.problem(), globI, /*timeIdx*/ 0);
             intQuants_[globI].updateEnergyQuantities_(simulator_.problem(), globI, /*timeIdx*/ 0);
         }
@@ -248,6 +362,9 @@ protected:
             }
         }
         else {
+            #ifdef _OPENMP
+            #pragma omp parallel for
+            #endif
             for (unsigned globI = 0; globI < numCells; ++globI) {
                 this->temperature_[globI] -= std::clamp(dx[globI][0], -this->maxTempChange_, this->maxTempChange_);
                 intQuants_[globI].updateTemperature_(simulator_.problem(), globI, /*timeIdx*/ 0);
@@ -266,6 +383,9 @@ protected:
         const IsNumericalAquiferCell isNumericalAquiferCell(simulator_.gridView().grid());
         Scalar sum_pv = 0.0;
         Scalar sum_pv_not_converged = 0.0;
+        #ifdef _OPENMP
+        #pragma omp parallel for
+        #endif
         for (const auto& elem : elements(simulator_.gridView(), Dune::Partitions::interior)) {
             unsigned globI = elemMapper.index(elem);
             const auto pvValue =  simulator_.problem().referencePorosity(globI, /*timeIdx=*/0)
@@ -313,9 +433,9 @@ protected:
     void computeStorageTerm(unsigned globI, LhsEval& storage)
     {
         const auto& intQuants = intQuants_[globI];
-        const auto& poro = decay<LhsEval>(intQuants.porosity());
+        const auto& poro = getValue(simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0).porosity());
         // accumulate the internal energy of the fluids
-        const auto& fs = intQuants.fluidState();
+        const auto& fs = intQuants.fluidStateTemp();
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
                 continue;
@@ -332,7 +452,7 @@ protected:
         const Scalar rockFraction = intQuants.rockFraction();
         const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
         storage += rockFraction*uRock;
-        storage*= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+        storage*= scalingFactor_;
     }
 
     template < class ResidualNBInfo>
@@ -340,8 +460,8 @@ protected:
                          const ResidualNBInfo& res_nbinfo,
                          Evaluation& flux)
     {
-        const IntensiveQuantities& intQuantsIn = intQuants_[globI];
-        const IntensiveQuantities& intQuantsEx = intQuants_[globJ];
+        const IntensiveQuantities& intQuantsIn = simulator_.model().intensiveQuantities(globI, /*timeIdx*/ 0);
+        const IntensiveQuantities& intQuantsEx = simulator_.model().intensiveQuantities(globJ, /*timeIdx*/ 0);
         RateVector tmp(0.0); //not used
         RateVector darcyFlux(0.0);
         LocalResidual::computeFlux(tmp, darcyFlux, globI, globJ, intQuantsIn, intQuantsEx,
@@ -355,12 +475,12 @@ protected:
                 FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
 
             bool inIsUp = darcyFlux[activeCompIdx] > 0;
-            const IntensiveQuantities& up = inIsUp ? intQuantsIn : intQuantsEx;
-            const auto& fs = up.fluidState();
+            const auto& up = inIsUp ? intQuants_[globI] : intQuants_[globJ];
+            const auto& fs = up.fluidStateTemp();
             if (inIsUp) {
                 flux += fs.enthalpy(phaseIdx)
                     * fs.density(phaseIdx)
-                    * darcyFlux[activeCompIdx];
+                    * getValue(darcyFlux[activeCompIdx]);
             }
             else {
                 flux += getValue(fs.enthalpy(phaseIdx))
@@ -368,7 +488,7 @@ protected:
                     * getValue(darcyFlux[activeCompIdx]);
             }
         }
-        flux *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+        flux *= scalingFactor_;
     }
 
     template < class ResidualNBInfo>
@@ -376,8 +496,8 @@ protected:
                              const ResidualNBInfo& res_nbinfo,
                              Evaluation& heatFlux)
     {
-        const IntensiveQuantities& intQuantsIn = intQuants_[globI];
-        const IntensiveQuantities& intQuantsEx = intQuants_[globJ];
+        const auto& intQuantsIn = intQuants_[globI];
+        const auto& intQuantsEx = intQuants_[globJ];
         const Scalar inAlpha = simulator_.problem().thermalHalfTransmissibility(globI, globJ);
         const Scalar outAlpha = simulator_.problem().thermalHalfTransmissibility(globJ, globI);
         short interiorDofIdx = 0; // NB
@@ -388,12 +508,12 @@ protected:
                                                         exteriorDofIdx,
                                                         intQuantsIn,
                                                         intQuantsEx,
-                                                        intQuantsIn.fluidState(),
-                                                        intQuantsEx.fluidState(),
+                                                        intQuantsIn.fluidStateTemp(),
+                                                        intQuantsEx.fluidStateTemp(),
                                                         inAlpha,
                                                         outAlpha,
                                                         res_nbinfo.faceArea);
-        heatFlux *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>()*res_nbinfo.faceArea;
+        heatFlux *= scalingFactor_*res_nbinfo.faceArea;
     }
 
     void assembleEquations()
@@ -411,7 +531,7 @@ protected:
             Evaluation storage = 0.0;
             computeStorageTerm(globI, storage);
             this->energyVector_[globI] += storefac * ( getValue(storage) - storage1_[globI] );
-            (*this->energyMatrix_)[globI][globI][0][0] += storefac * storage.derivative(Indices::temperatureIdx);
+            (*this->energyMatrix_)[globI][globI][0][0] += storefac * storage.derivative(temperatureIdx);
         }
 
         const auto& neighborInfo = simulator_.model().linearizer().getNeighborInfo();
@@ -428,15 +548,15 @@ protected:
                 Evaluation flux = 0.0;
                 computeFluxTerm(globI, globJ, nbInfo.res_nbinfo, flux);
                 this->energyVector_[globI] += getValue(flux);
-                (*this->energyMatrix_)[globI][globI][0][0] += flux.derivative(Indices::temperatureIdx);
-                (*this->energyMatrix_)[globJ][globI][0][0] -= flux.derivative(Indices::temperatureIdx);
+                (*this->energyMatrix_)[globI][globI][0][0] += flux.derivative(temperatureIdx);
+                (*this->energyMatrix_)[globJ][globI][0][0] -= flux.derivative(temperatureIdx);
 
                 // compute conductive flux
                 Evaluation heatFlux = 0.0;
                 computeHeatFluxTerm(globI, globJ, nbInfo.res_nbinfo, heatFlux);
                 this->energyVector_[globI] += getValue(heatFlux);
-                (*this->energyMatrix_)[globI][globI][0][0] += heatFlux.derivative(Indices::temperatureIdx);
-                (*this->energyMatrix_)[globJ][globI][0][0] -= heatFlux.derivative(Indices::temperatureIdx);
+                (*this->energyMatrix_)[globI][globI][0][0] += heatFlux.derivative(temperatureIdx);
+                (*this->energyMatrix_)[globJ][globI][0][0] -= heatFlux.derivative(temperatureIdx);
             }
         }
 
@@ -450,14 +570,17 @@ protected:
 
         bool enableDriftCompensation = Parameters::Get<Parameters::EnableDriftCompensationTemp>();
         if (enableDriftCompensation) {
+            #ifdef _OPENMP
+            #pragma omp parallel for
+            #endif
             for (unsigned globalDofIdx = 0; globalDofIdx < numCells; ++globalDofIdx) {
                 auto dofDriftRate = problem.drift()[globalDofIdx]/dt;
-                const auto& fs = intQuants_[globalDofIdx].fluidState();
+                const auto& fs = intQuants_[globalDofIdx].fluidStateTemp();
                 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
                    const unsigned activeCompIdx =
                         FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
                    auto drift_hrate = dofDriftRate[activeCompIdx]*getValue(fs.enthalpy(phaseIdx)) * getValue(fs.density(phaseIdx)) /  getValue(fs.invB(phaseIdx));
-                   this->energyVector_[globalDofIdx] -= drift_hrate*getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+                   this->energyVector_[globalDofIdx] -= drift_hrate*scalingFactor_;
                 }
             }
         }
@@ -484,7 +607,7 @@ protected:
         this->energy_rates_[well_index] = 0.0;
         for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
             const auto globI = ws.perf_data.cell_index[i];
-            auto fs = intQuants_[globI].fluidState(); //copy to make it possible to change the temp in the injector
+            auto fs = intQuants_[globI].fluidStateTemp(); //copy to make it possible to change the temp in the injector
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
                 if (!FluidSystem::phaseIsActive(phaseIdx)) {
                     continue;
@@ -513,18 +636,19 @@ protected:
                     rate *= fs.enthalpy(phaseIdx) * getValue(fs.density(phaseIdx)) / getValue(fs.invB(phaseIdx));
                 }
                 this->energy_rates_[well_index] += getValue(rate);
-                rate *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+                rate *= scalingFactor_;
                 this->energyVector_[globI] -= getValue(rate);
-                (*this->energyMatrix_)[globI][globI][0][0] -= rate.derivative(Indices::temperatureIdx);
+                (*this->energyMatrix_)[globI][globI][0][0] -= rate.derivative(temperatureIdx);
             }
         }
     }
 
     const Simulator& simulator_;
     EnergyVector storage1_;
-    std::vector<IntensiveQuantities> intQuants_;
+    std::vector<IntensiveQuantitiesTemp> intQuants_;
     std::vector<int> overlapRows_;
     std::vector<int> interiorRows_;
+    Scalar scalingFactor_{1.0};
 };
 
 // need for the old linearizer

opm/simulators/utils/InstantiationIndicesMacros.hpp

@@ -63,30 +63,29 @@
     template class CLASS<FS<T>,__VA_ARGS__>;
 
 #define INSTANTIATE_TYPE_INDICES(CLASS, T)                                                  \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,false,0u,1u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,false,0u,1u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,false,0u,1u,5u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,0u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,1u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,2u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,true,0u,2u,0u>)  \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,false,0u,2u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,false,0u,2u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,false,0u,1u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,true,0u,0u,0u>)  \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,false,0u,0u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,true,0u,0u,0u>)  \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<1u,0u,0u,0u,false,false,false,0u,0u,0u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,false,0u,0u,2u>) \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,false,false,0u,0u>)            \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,true,false,false,0u,0u>)            \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,true,false,0u,0u>)             \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,false,true,0u,0u>)             \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,false,false,false,0u,0u>)            \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,1u,0u,0u,false,false,false,0u,0u>)            \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,1u,0u,false,false,false,0u,0u>)            \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,false,false,0u,0u>)            \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,false,true,0u,0u>)             \
-    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,false,true,false,0u,0u>)
+    INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,5u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u,0u>)  \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,1u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,0u,0u>)  \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,false,0u,0u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,1u,false,true,0u,0u,0u>)  \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<1u,0u,0u,0u,false,false,0u,0u,0u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u,2u>) \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,false,0u,0u>)            \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,true,false,0u,0u>)             \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,0u,false,true,0u,0u>)             \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,false,false,0u,0u>)            \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,1u,0u,0u,false,false,0u,0u>)            \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,1u,0u,false,false,0u,0u>)            \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,false,0u,0u>)            \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<0u,0u,0u,1u,false,true,0u,0u>)             \
+    INSTANTIATE_CLASS(CLASS,T,BlackOilVariableAndEquationIndices<1u,0u,0u,0u,true,false,0u,0u>)
 
 #endif

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -783,7 +783,13 @@ namespace Opm {
         this->commitWGState();
 
         //reporting output temperatures
-        this->computeWellTemperature();
+        const auto& eclState = simulator_.vanguard().eclState();
+        bool isTemp = eclState.getSimulationConfig().isTemp();
+        bool isThermal = eclState.getSimulationConfig().isThermal();
+        if (isTemp || isThermal) {
+            this->computeWellTemperature();
+        }
+
     }
 
 
@@ -2173,8 +2179,7 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     computeWellTemperature()
     {
-        if constexpr (energyModuleType_ == EnergyModules::FullyImplicitThermal ||
-                      energyModuleType_ == EnergyModules::SequentialImplicitThermal) {
+        if constexpr (energyModuleType_ == EnergyModules::FullyImplicitThermal) {
             int np = this->numPhases();
             Scalar cellInternalEnergy;
             Scalar cellBinv;

opm/simulators/wells/WellInterface.hpp

@@ -129,8 +129,8 @@ public:
     // For the conversion between the surface volume rate and reservoir voidage rate
     using FluidState = BlackOilFluidState<Eval,
                                           FluidSystem,
-                                          energyModuleType == EnergyModules::ConstantTemperature,
-                                          (energyModuleType == EnergyModules::FullyImplicitThermal || energyModuleType == EnergyModules::SequentialImplicitThermal),
+                                          energyModuleType != EnergyModules::NoTemperature,
+                                          energyModuleType == EnergyModules::FullyImplicitThermal,
                                           Indices::compositionSwitchIdx >= 0,
                                           has_watVapor,
                                           has_brine,

python/simulators/PyGasWaterSimulator.cpp

@@ -46,7 +46,6 @@ namespace Opm {
                                                      getPropValue<TypeTag, Properties::EnableExtbo>(),
                                                      getPropValue<TypeTag, Properties::EnablePolymer>(),
                                                      getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal,
-                                                     getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal,
                                                      getPropValue<TypeTag, Properties::EnableFoam>(),
                                                      getPropValue<TypeTag, Properties::EnableBrine>(),
                                                      /*PVOffset=*/0,

python/simulators/PyOnePhaseSimulator.cpp

@@ -45,7 +45,6 @@ namespace Opm::Properties {
                                              getPropValue<TypeTag, Properties::EnableExtbo>(),
                                              getPropValue<TypeTag, Properties::EnablePolymer>(),
                                              getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal,
-                                             getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::SequentialImplicitThermal,
                                              getPropValue<TypeTag, Properties::EnableFoam>(),
                                              getPropValue<TypeTag, Properties::EnableBrine>(),
                                              /*PVOffset=*/0,