Merge pull request #822 from hnil/tpfa_temp

TPFA version for energy

opm/models/blackoil/blackoilenergymodules.hh

@@ -31,6 +31,7 @@
 #include "blackoilproperties.hh"
 #include <opm/models/io/vtkblackoilenergymodule.hh>
 #include <opm/models/common/quantitycallbacks.hh>
+#include <opm/models/discretization/common/linearizationtype.hh>
 
 #include <opm/material/common/Tabulated1DFunction.hpp>
 
@@ -53,7 +54,6 @@ class BlackOilEnergyModule
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
     using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
-    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
     using Model = GetPropType<TypeTag, Properties::Model>;
     using Simulator = GetPropType<TypeTag, Properties::Simulator>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
@@ -70,6 +70,7 @@ class BlackOilEnergyModule
     static constexpr unsigned numPhases = FluidSystem::numPhases;
 
 public:
+    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
     /*!
      * \brief Register all run-time parameters for the black-oil energy module.
      */
@@ -191,6 +192,30 @@ public:
         }
     }
 
+    static void addHeatFlux(RateVector& flux,
+                            const Evaluation& heatFlux)
+    {
+        if constexpr (enableEnergy) {
+            // diffusive energy flux
+            flux[contiEnergyEqIdx] += heatFlux;
+            flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+        }
+    }
+
+
+
+    template <class UpEval, class Eval, class FluidState>
+    static void addPhaseEnthalpyFluxes_(RateVector& flux,
+                                        unsigned phaseIdx,
+                                        const Eval& volumeFlux,
+                                        const FluidState& upFs)
+    {
+        flux[contiEnergyEqIdx] +=
+            decay<UpEval>(upFs.enthalpy(phaseIdx))
+            * decay<UpEval>(upFs.density(phaseIdx))
+            * volumeFlux;
+    }
+
     template <class UpstreamEval>
     static void addPhaseEnthalpyFlux_(RateVector& flux,
                                       unsigned phaseIdx,
@@ -202,12 +227,11 @@ public:
         unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
         const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
         const auto& fs = up.fluidState();
-
         const auto& volFlux = extQuants.volumeFlux(phaseIdx);
-        flux[contiEnergyEqIdx] +=
+        addPhaseEnthalpyFluxes_<UpstreamEval>(flux,
-            decay<UpstreamEval>(fs.enthalpy(phaseIdx))
+                                              phaseIdx,
-            * decay<UpstreamEval>(fs.density(phaseIdx))
+                                              volFlux,
-            * volFlux;
+                                              fs);
     }
 
     static void addToEnthalpyRate(RateVector& flux,
@@ -317,6 +341,7 @@ class BlackOilEnergyIntensiveQuantities
     using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
 
@@ -341,6 +366,20 @@ public:
         fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, elemCtx.linearizationType()));
     }
 
+    /*!
+     * \brief Update the temperature of the intensive quantity's fluid state
+     *
+     */
+    void updateTemperature_([[maybe_unused]] const Problem& problem,
+                            const PrimaryVariables& priVars,
+                            [[maybe_unused]] unsigned globalDofIdx,
+                            const unsigned timeIdx,
+                            const LinearizationType& lintype)
+    {
+        auto& fs = asImp_().fluidState_;
+        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, lintype));
+    }
+
     /*!
      * \brief Compute the intensive quantities needed to handle energy conservation
      *
@@ -400,12 +439,13 @@ template <class TypeTag>
 class BlackOilEnergyIntensiveQuantities<TypeTag, false>
 {
     using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
-
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
     static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>();
 
 public:
@@ -417,7 +457,24 @@ public:
             // even if energy is conserved, the temperature can vary over the spatial
             // domain if the EnableTemperature property is set to true
             auto& fs = asImp_().fluidState_;
-            Scalar T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
+            const Scalar T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
+            fs.setTemperature(T);
+        }
+    }
+
+    template<class Problem>
+    void updateTemperature_([[maybe_unused]] const Problem& problem,
+                            [[maybe_unused]] const PrimaryVariables& priVars,
+                            [[maybe_unused]] unsigned globalDofIdx,
+                            [[maybe_unused]] unsigned timeIdx,
+                            [[maybe_unused]] const LinearizationType& lintype
+        )
+    {
+        if constexpr (enableTemperature) {
+            auto& fs = asImp_().fluidState_;
+            // even if energy is conserved, the temperature can vary over the spatial
+            // domain if the EnableTemperature property is set to true
+            const Scalar T = problem.temperature(globalDofIdx, timeIdx);
             fs.setTemperature(T);
         }
     }
@@ -469,29 +526,26 @@ class BlackOilEnergyExtensiveQuantities
     static const int dimWorld = GridView::dimensionworld;
     using DimVector = Dune::FieldVector<Scalar, dimWorld>;
     using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
-
 public:
-    void updateEnergy(const ElementContext& elemCtx,
+    template<class FluidState>
-                      unsigned scvfIdx,
+    static void updateEnergy(Evaluation& energyFlux,
-                      unsigned timeIdx)
+                             const unsigned& focusDofIndex,
+                             const unsigned& inIdx,
+                             const unsigned& exIdx,
+                             const IntensiveQuantities& inIq,
+                             const IntensiveQuantities& exIq,
+                             const FluidState& inFs,
+                             const FluidState& exFs,
+                             const Scalar& inAlpha,
+                             const Scalar& outAlpha,
+                             const Scalar& faceArea)
     {
-        const auto& stencil = elemCtx.stencil(timeIdx);
-        const auto& scvf = stencil.interiorFace(scvfIdx);
-
-        Scalar faceArea = scvf.area();
-        unsigned inIdx = scvf.interiorIndex();
-        unsigned exIdx = scvf.exteriorIndex();
-        const auto& inIq = elemCtx.intensiveQuantities(inIdx, timeIdx);
-        const auto& exIq = elemCtx.intensiveQuantities(exIdx, timeIdx);
-        const auto& inFs = inIq.fluidState();
-        const auto& exFs = exIq.fluidState();
-
         Evaluation deltaT;
-        if (elemCtx.focusDofIndex() == inIdx)
+        if (focusDofIndex == inIdx)
             deltaT =
                 decay<Scalar>(exFs.temperature(/*phaseIdx=*/0))
                 - inFs.temperature(/*phaseIdx=*/0);
-        else if (elemCtx.focusDofIndex() == exIdx)
+        else if (focusDofIndex == exIdx)
             deltaT =
                 exFs.temperature(/*phaseIdx=*/0)
                 - decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
@@ -501,28 +555,23 @@ public:
                 - decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
 
         Evaluation inLambda;
-        if (elemCtx.focusDofIndex() == inIdx)
+        if (focusDofIndex == inIdx)
             inLambda = inIq.totalThermalConductivity();
         else
             inLambda = decay<Scalar>(inIq.totalThermalConductivity());
 
         Evaluation exLambda;
-        if (elemCtx.focusDofIndex() == exIdx)
+        if (focusDofIndex == exIdx)
             exLambda = exIq.totalThermalConductivity();
         else
             exLambda = decay<Scalar>(exIq.totalThermalConductivity());
 
-        auto distVec = elemCtx.pos(exIdx, timeIdx);
-        distVec -= elemCtx.pos(inIdx, timeIdx);
-
         Evaluation H;
         if (inLambda > 0.0 && exLambda > 0.0) {
             // compute the "thermal transmissibility". In contrast to the normal
             // transmissibility this cannot be done as a preprocessing step because the
-            // average thermal thermal conductivity is analogous to the permeability but
+            // average thermal conductivity is analogous to the permeability but
             // depends on the solution.
-            Scalar inAlpha = elemCtx.problem().thermalHalfTransmissibilityIn(elemCtx, scvfIdx, timeIdx);
-            Scalar outAlpha = elemCtx.problem().thermalHalfTransmissibilityOut(elemCtx, scvfIdx, timeIdx);
             const Evaluation& inH = inLambda*inAlpha;
             const Evaluation& exH = exLambda*outAlpha;
             H = 1.0/(1.0/inH + 1.0/exH);
@@ -530,7 +579,36 @@ public:
         else
             H = 0.0;
 
-        energyFlux_ = deltaT * (-H/faceArea);
+        energyFlux = deltaT * (-H/faceArea);
+    }
+
+    void updateEnergy(const ElementContext& elemCtx,
+                      unsigned scvfIdx,
+                      unsigned timeIdx)
+    {
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& scvf = stencil.interiorFace(scvfIdx);
+
+        const Scalar faceArea = scvf.area();
+        const unsigned inIdx = scvf.interiorIndex();
+        const unsigned exIdx = scvf.exteriorIndex();
+        const auto& inIq = elemCtx.intensiveQuantities(inIdx, timeIdx);
+        const auto& exIq = elemCtx.intensiveQuantities(exIdx, timeIdx);
+        const auto& inFs = inIq.fluidState();
+        const auto& exFs = exIq.fluidState();
+        const Scalar inAlpha = elemCtx.problem().thermalHalfTransmissibilityIn(elemCtx, scvfIdx, timeIdx);
+        const Scalar outAlpha = elemCtx.problem().thermalHalfTransmissibilityOut(elemCtx, scvfIdx, timeIdx);
+        updateEnergy(energyFlux_,
+                     elemCtx.focusDofIndex(),
+                     inIdx,
+                     exIdx,
+                     inIq,
+                     exIq,
+                     inFs,
+                     exFs,
+                     inAlpha,
+                     outAlpha,
+                     faceArea);
     }
 
     template <class Context, class BoundaryFluidState>
@@ -544,10 +622,22 @@ public:
 
         unsigned inIdx = scvf.interiorIndex();
         const auto& inIq = ctx.intensiveQuantities(inIdx, timeIdx);
-        const auto& inFs = inIq.fluidState();
+        const auto& focusDofIdx = ctx.focusDofIndex();
+        const Scalar alpha = ctx.problem().thermalHalfTransmissibilityBoundary(ctx, scvfIdx);
+        updateEnergyBoundary(energyFlux_, inIq, focusDofIdx, inIdx, alpha, boundaryFs);
+    }
 
+    template <class BoundaryFluidState>
+    static void updateEnergyBoundary(Evaluation& energyFlux,
+                                     const IntensiveQuantities& inIq,
+                                     unsigned focusDofIndex,
+                                     unsigned inIdx,
+                                     Scalar alpha,
+                                     const BoundaryFluidState& boundaryFs)
+    {
+        const auto& inFs = inIq.fluidState();
         Evaluation deltaT;
-        if (ctx.focusDofIndex() == inIdx)
+        if (focusDofIndex == inIdx)
             deltaT =
                 boundaryFs.temperature(/*phaseIdx=*/0)
                 - inFs.temperature(/*phaseIdx=*/0);
@@ -557,24 +647,21 @@ public:
                 - decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
 
         Evaluation lambda;
-        if (ctx.focusDofIndex() == inIdx)
+        if (focusDofIndex == inIdx)
             lambda = inIq.totalThermalConductivity();
         else
             lambda = decay<Scalar>(inIq.totalThermalConductivity());
 
-        auto distVec = scvf.integrationPos();
-        distVec -= ctx.pos(inIdx, timeIdx);
 
         if (lambda > 0.0) {
             // compute the "thermal transmissibility". In contrast to the normal
             // transmissibility this cannot be done as a preprocessing step because the
             // average thermal conductivity is analogous to the permeability but depends
             // on the solution.
-            Scalar alpha = ctx.problem().thermalHalfTransmissibilityBoundary(ctx, scvfIdx);
+            energyFlux = deltaT*lambda*(-alpha);
-            energyFlux_ = deltaT*lambda*(-alpha);
         }
         else
-            energyFlux_ = 0.0;
+            energyFlux = 0.0;
     }
 
     const Evaluation& energyFlux()  const
@@ -592,8 +679,23 @@ class BlackOilEnergyExtensiveQuantities<TypeTag, false>
 {
     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>
+    static void updateEnergy(Evaluation& energyFlux,
+                             const unsigned& focusDofIndex,
+                             const unsigned& inIdx,
+                             const unsigned& exIdx,
+                             const IntensiveQuantities& inIq,
+                             const IntensiveQuantities& exIq,
+                             const FluidState& inFs,
+                             const FluidState& exFs,
+                             const Scalar& inAlpha,
+                             const Scalar& outAlpha,
+                             const Scalar& faceArea)
+    {};
+
     void updateEnergy(const ElementContext&,
                       unsigned,
                       unsigned)
@@ -606,6 +708,16 @@ public:
                               const BoundaryFluidState&)
     {}
 
+    template <class BoundaryFluidState>
+    static void updateEnergyBoundary(Evaluation& heatFlux,
+                                     const IntensiveQuantities& inIq,
+                                     unsigned focusDofIndex,
+                                     unsigned inIdx,
+                                     unsigned timeIdx,
+                                     Scalar alpha,
+                                     const BoundaryFluidState& boundaryFs){
+    }
+
     const Evaluation& energyFlux()  const
     { throw std::logic_error("Requested the energy flux, but energy is not conserved"); }
 };

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -106,6 +106,19 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     using Toolbox = MathToolbox<Evaluation>;
 
 public:
+
+    struct ResidualNBInfo
+    {
+        double trans;
+        double faceArea;
+        double thpres;
+        double dZg;
+        FaceDir::DirEnum faceDirection;
+        double Vin;
+        double Vex;
+        double inAlpha;
+        double outAlpha;
+    };
     /*!
      * \copydoc FvBaseLocalResidual::computeStorage
      */
@@ -205,53 +218,23 @@ public:
      */
     static void computeFlux(RateVector& flux,
                             RateVector& darcy,
-                            const Problem& problem,
                             const unsigned globalIndexIn,
                             const unsigned globalIndexEx,
                             const IntensiveQuantities& intQuantsIn,
                             const IntensiveQuantities& intQuantsEx,
-                            const Scalar trans,
+                            const ResidualNBInfo& nbInfo)
-                            const Scalar faceArea,
-                            const FaceDir::DirEnum facedir)
     {
         OPM_TIMEBLOCK_LOCAL(computeFlux);
         flux = 0.0;
         darcy = 0.0;
-        Scalar Vin = problem.model().dofTotalVolume(globalIndexIn);
-        Scalar Vex = problem.model().dofTotalVolume(globalIndexEx);
-
-        Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx);
-
-        // estimate the gravity correction: for performance reasons we use a simplified
-        // approach for this flux module that assumes that gravity is constant and always
-        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
-        Scalar g = problem.gravity()[dimWorld - 1];
-
-        // this is quite hacky because the dune grid interface does not provide a
-        // cellCenterDepth() method (so we ask the problem to provide it). The "good"
-        // solution would be to take the Z coordinate of the element centroids, but since
-        // ECL seems to like to be inconsistent on that front, it needs to be done like
-        // here...
-        Scalar zIn = problem.dofCenterDepth(globalIndexIn);
-        Scalar zEx = problem.dofCenterDepth(globalIndexEx);
-
-        // the distances from the DOF's depths. (i.e., the additional depth of the
-        // exterior DOF)
-        Scalar distZ = zIn - zEx; // NB could be precalculated
 
         calculateFluxes_(flux,
                          darcy,
                          intQuantsIn,
                          intQuantsEx,
-                         Vin,
-                         Vex,
                          globalIndexIn,
                          globalIndexEx,
-                         distZ * g,
+                         nbInfo);
-                         thpres,
-                         trans,
-                         faceArea,
-                         facedir);
     }
 
     // This function demonstrates compatibility with the ElementContext-based interface.
@@ -262,7 +245,7 @@ public:
                             unsigned scvfIdx,
                             unsigned timeIdx)
     {
-        OPM_TIMEBLOCK_LOCAL(computeFlux);              
+        OPM_TIMEBLOCK_LOCAL(computeFlux);
         assert(timeIdx == 0);
 
         flux = 0.0;
@@ -303,43 +286,45 @@ public:
         // solution would be to take the Z coordinate of the element centroids, but since
         // ECL seems to like to be inconsistent on that front, it needs to be done like
         // here...
-        Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
+        const Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
-        Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
+        const Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
-
         // the distances from the DOF's depths. (i.e., the additional depth of the
         // exterior DOF)
-        Scalar distZ = zIn - zEx;
+        const Scalar distZ = zIn - zEx;
+        // for thermal harmonic mean of half trans
+        const Scalar inAlpha = problem.thermalHalfTransmissibility(globalIndexIn, globalIndexEx);
+        const Scalar outAlpha = problem.thermalHalfTransmissibility(globalIndexEx, globalIndexIn);
+
+        const ResidualNBInfo res_nbinfo {trans, faceArea, thpres, distZ * g, facedir, Vin, Vex, inAlpha, outAlpha};
 
         calculateFluxes_(flux,
                          darcy,
                          intQuantsIn,
                          intQuantsEx,
-                         Vin,
-                         Vex,
                          globalIndexIn,
                          globalIndexEx,
-                         distZ * g,
+                         res_nbinfo);
-                         thpres,
-                         trans,
-                         faceArea,
-                         facedir);
     }
 
     static void calculateFluxes_(RateVector& flux,
                                  RateVector& darcy,
                                  const IntensiveQuantities& intQuantsIn,
                                  const IntensiveQuantities& intQuantsEx,
-                                 const Scalar& Vin,
-                                 const Scalar& Vex,
                                  const unsigned& globalIndexIn,
                                  const unsigned& globalIndexEx,
-                                 const Scalar& distZg,
+                                 const ResidualNBInfo& nbInfo)
-                                 const Scalar& thpres,
-                                 const Scalar& trans,
-                                 const Scalar& faceArea,
-                                 const FaceDir::DirEnum facedir)
     {
         OPM_TIMEBLOCK_LOCAL(calculateFluxes);
+        const Scalar Vin = nbInfo.Vin;
+        const Scalar Vex = nbInfo.Vex;
+        const Scalar distZg = nbInfo.dZg;
+        const Scalar thpres = nbInfo.thpres;
+        const Scalar trans = nbInfo.trans;
+        const Scalar faceArea = nbInfo.faceArea;
+        const FaceDir::DirEnum facedir = nbInfo.faceDirection;
+        const Scalar inAlpha = nbInfo.inAlpha;
+        const Scalar outAlpha = nbInfo.outAlpha;
+
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx))
                 continue;
@@ -358,7 +343,7 @@ public:
                                                              intQuantsEx,
                                                              phaseIdx, // input
                                                              interiorDofIdx, // input
-                                                             exteriorDofIdx, // intput
+                                                             exteriorDofIdx, // input
                                                              Vin,
                                                              Vex,
                                                              globalIndexIn,
@@ -382,7 +367,7 @@ public:
                        (Toolbox::value(up.mobility(phaseIdx, facedir)) * Toolbox::value(transMult) * (-trans / faceArea));
             }
             unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
-            darcy[conti0EqIdx + activeCompIdx] = darcyFlux.value() * faceArea; // For the FLORES fluxes
+            darcy[conti0EqIdx + activeCompIdx] = darcyFlux.value() * faceArea; // NB! For the FLORES fluxes without derivatives
 
             unsigned pvtRegionIdx = up.pvtRegionIndex();
             // if (upIdx == globalFocusDofIdx){
@@ -392,12 +377,22 @@ public:
                 const auto& surfaceVolumeFlux = invB * darcyFlux;
                 evalPhaseFluxes_<Evaluation, Evaluation, FluidState>(
                     flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
+                if constexpr (enableEnergy) {
+                    EnergyModule::template addPhaseEnthalpyFluxes_<Evaluation, Evaluation, FluidState>(
+                        flux, phaseIdx, darcyFlux, up.fluidState());
+                }
             } else {
                 const auto& invB = getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), phaseIdx, pvtRegionIdx);
                 const auto& surfaceVolumeFlux = invB * darcyFlux;
                 evalPhaseFluxes_<Scalar, Evaluation, FluidState>(
                     flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
+                if constexpr (enableEnergy) {
+                    EnergyModule::template
+                        addPhaseEnthalpyFluxes_<Scalar, Evaluation, FluidState>
+                        (flux,phaseIdx,darcyFlux, up.fluidState());
+                }
             }
+
         }
 
         // deal with solvents (if present)
@@ -413,7 +408,27 @@ public:
         // PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
         // deal with energy (if present)
-        static_assert(!enableEnergy, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        if constexpr(enableEnergy){
+            Evaluation heatFlux;
+            {
+                short interiorDofIdx = 0; // NB
+                short exteriorDofIdx = 1; // NB
+
+                EnergyModule::ExtensiveQuantities::template updateEnergy(heatFlux,
+                                                                         interiorDofIdx, // focusDofIndex,
+                                                                         interiorDofIdx,
+                                                                         exteriorDofIdx,
+                                                                         intQuantsIn,
+                                                                         intQuantsEx,
+                                                                         intQuantsIn.fluidState(),
+                                                                         intQuantsEx.fluidState(),
+                                                                         inAlpha,
+                                                                         outAlpha,
+                                                                         faceArea);
+            }
+            EnergyModule::addHeatFlux(flux, heatFlux);
+        }
+        // NB need to be tha last energy call since it does scaling
         // EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
         // deal with foam (if present)
@@ -494,37 +509,63 @@ public:
             const unsigned pvtRegionIdx = insideIntQuants.pvtRegionIndex();
 
             RateVector tmp;
-
+            const auto& darcyFlux = volumeFlux[phaseIdx];
             // mass conservation
             if (pBoundary < pInside) {
                 // outflux
                 const auto& invB = getInvB_<FluidSystem, FluidState, Evaluation>(insideIntQuants.fluidState(), phaseIdx, pvtRegionIdx);
-                Evaluation surfaceVolumeFlux = invB * volumeFlux[phaseIdx];
+                Evaluation surfaceVolumeFlux = invB * darcyFlux;
                 evalPhaseFluxes_<Evaluation>(tmp,
                                              phaseIdx,
                                              insideIntQuants.pvtRegionIndex(),
                                              surfaceVolumeFlux,
                                              insideIntQuants.fluidState());
+                if constexpr (enableEnergy) {
+                    EnergyModule::template
+                        addPhaseEnthalpyFluxes_<Evaluation, Evaluation, FluidState>
+                        (tmp, phaseIdx, darcyFlux, insideIntQuants.fluidState());
+                }
             } else if (pBoundary > pInside) {
                 // influx
                 using ScalarFluidState = decltype(bdyInfo.exFluidState);
                 const auto& invB = getInvB_<FluidSystem, ScalarFluidState, Scalar>(bdyInfo.exFluidState, phaseIdx, pvtRegionIdx);
-                Evaluation surfaceVolumeFlux = invB * volumeFlux[phaseIdx];
+                Evaluation surfaceVolumeFlux = invB * darcyFlux;
                 evalPhaseFluxes_<Scalar>(tmp,
                                          phaseIdx,
                                          insideIntQuants.pvtRegionIndex(),
                                          surfaceVolumeFlux,
                                          bdyInfo.exFluidState);
+                if constexpr (enableEnergy) {
+                    EnergyModule::template
+                        addPhaseEnthalpyFluxes_<Scalar, Evaluation, ScalarFluidState>
+                        (tmp,
+                         phaseIdx,
+                         darcyFlux,
+                         bdyInfo.exFluidState);
+                }
             }
 
             for (unsigned i = 0; i < tmp.size(); ++i) {
                 bdyFlux[i] += tmp[i];
             }
-
-            static_assert(!enableEnergy, "Relevant treatment of boundary conditions must be implemented before enabling.");
-            // Add energy flux treatment per phase here.
         }
 
+        // conductive heat flux from boundary
+        Evaluation heatFlux;
+        if constexpr(enableEnergy){
+            // avoid overload of functions with same numeber of elements in eclproblem
+            Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
+            unsigned inIdx = 0;//dummy
+            // always calculated with derivatives of this cell
+            EnergyModule::ExtensiveQuantities::template updateEnergyBoundary(heatFlux,
+                                                                             insideIntQuants,
+                                                                             /*focusDofIndex*/ inIdx,
+                                                                             inIdx,
+                                                                             alpha,
+                                                                             bdyInfo.exFluidState);
+        }
+        EnergyModule::addHeatFlux(bdyFlux, heatFlux);
+
         static_assert(!enableSolvent, "Relevant treatment of boundary conditions must be implemented before enabling.");
         static_assert(!enablePolymer, "Relevant treatment of boundary conditions must be implemented before enabling.");
         static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling.");
@@ -532,9 +573,6 @@ public:
         // make sure that the right mass conservation quantities are used
         adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex());
 
-        // heat conduction
-        static_assert(!enableEnergy, "Relevant treatment of boundary conditions must be implemented before enabling.");
-
 #ifndef NDEBUG
         for (unsigned i = 0; i < numEq; ++i) {
             Valgrind::CheckDefined(bdyFlux[i]);
@@ -596,8 +634,8 @@ public:
         MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
 
         // scale the source term of the energy equation
-        if (enableEnergy)
+        if constexpr(enableEnergy)
-            source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+                        source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
     template <class UpEval, class FluidState>

opm/models/blackoil/blackoilmodel.hh

@@ -139,11 +139,9 @@ struct Indices<TypeTag, TTag::BlackOilModel>
 template<class TypeTag>
 struct FluidSystem<TypeTag, TTag::BlackOilModel>
 {
-private:
+public:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
-
-public:
     using type = BlackOilFluidSystem<Scalar>;
 };
 

opm/models/discretization/common/tpfalinearizer.hh

@@ -111,7 +111,7 @@ class TpfaLinearizer
     using ADVectorBlock = GetPropType<TypeTag, Properties::RateVector>;
 
     static const bool linearizeNonLocalElements = getPropValue<TypeTag, Properties::LinearizeNonLocalElements>();
-
+    static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
     // copying the linearizer is not a good idea
     TpfaLinearizer(const TpfaLinearizer&);
 //! \endcond
@@ -316,7 +316,7 @@ public:
     const auto& getFlowsInfo() const{
 
         return flowsInfo_;
-    }   
+    }
 
     /*!
      * \brief Return constant reference to the floresInfo.
@@ -419,7 +419,7 @@ private:
         // freedom of each primary degree of freedom
         using NeighborSet = std::set< unsigned >;
         std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
-
+        const Scalar gravity = problem_().gravity()[dimWorld - 1];
         unsigned numCells = model.numTotalDof();
         neighborInfo_.reserve(numCells, 6 * numCells);
         std::vector<NeighborInfo> loc_nbinfo;
@@ -436,15 +436,27 @@ private:
                     unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
                     sparsityPattern[myIdx].insert(neighborIdx);
                     if (dofIdx > 0) {
-                        const double trans = problem_().transmissibility(myIdx, neighborIdx);
+                        const Scalar trans = problem_().transmissibility(myIdx, neighborIdx);
                         const auto scvfIdx = dofIdx - 1;
                         const auto& scvf = stencil.interiorFace(scvfIdx);
-                        const double area = scvf.area();
+                        const Scalar area = scvf.area();
+                        const Scalar Vin = problem_().model().dofTotalVolume(myIdx);
+                        const Scalar Vex = problem_().model().dofTotalVolume(neighborIdx);
+                        const Scalar zIn = problem_().dofCenterDepth(myIdx);
+                        const Scalar zEx = problem_().dofCenterDepth(neighborIdx);
+                        const Scalar dZg = (zIn - zEx)*gravity;
+                        const Scalar thpres = problem_().thresholdPressure(myIdx, neighborIdx);
+                        Scalar inAlpha {0.};
+                        Scalar outAlpha {0.};
                         FaceDirection dirId = FaceDirection::Unknown;
+                        if constexpr(enableEnergy){
+                            inAlpha = problem_().thermalHalfTransmissibility(myIdx, neighborIdx);
+                            outAlpha = problem_().thermalHalfTransmissibility(neighborIdx, myIdx);
+                        }
                         if (materialLawManager->hasDirectionalRelperms()) {
                             dirId = scvf.faceDirFromDirId();
                         }
-                        loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area, dirId, nullptr};
+                        loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, {trans, area, thpres, dZg, dirId, Vin, Vex, inAlpha, outAlpha}, nullptr};
                     }
                 }
                 neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
@@ -551,7 +563,7 @@ private:
                         const int cartMyIdx = simulator_().vanguard().cartesianIndex(myIdx);
                         const int cartNeighborIdx = simulator_().vanguard().cartesianIndex(neighborIdx);
                         const auto& range = nncIndices.equal_range(cartMyIdx);
-                        for (auto it = range.first; it != range.second; ++it) { 
+                        for (auto it = range.first; it != range.second; ++it) {
                             if (it->second.first == cartNeighborIdx){
                                 // -1 gives problem since is used for the nncInput from the deck
                                 faceId = -2;
@@ -628,7 +640,7 @@ private:
 
             // Flux term.
             {
-            OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell);    
+            OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachCell);
             short loc = 0;
             for (const auto& nbInfo : nbInfos) {
                 OPM_TIMEBLOCK_LOCAL(fluxCalculationForEachFace);
@@ -639,10 +651,8 @@ private:
                 adres = 0.0;
                 darcyFlux = 0.0;
                 const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
-                LocalResidual::computeFlux(
+                LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo);
-                       adres, darcyFlux, problem_(), globI, globJ, intQuantsIn, intQuantsEx,
+                adres *= nbInfo.res_nbinfo.faceArea;
-                           nbInfo.trans, nbInfo.faceArea, nbInfo.faceDirection);
-                adres *= nbInfo.faceArea;
                 if (enableFlows) {
                     for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
                         flowsInfo_[globI][loc].flow[phaseIdx] = adres[phaseIdx].value();
@@ -773,7 +783,7 @@ private:
             auto nbInfos = neighborInfo_[globI]; // nbInfos will be a SparseTable<...>::mutable_iterator_range.
             for (auto& nbInfo : nbInfos) {
                 unsigned globJ = nbInfo.neighbor;
-                nbInfo.trans = problem_().transmissibility(globI, globJ);
+                nbInfo.res_nbinfo.trans = problem_().transmissibility(globI, globJ);
             }
         }
     }
@@ -789,12 +799,11 @@ private:
 
     LinearizationType linearizationType_;
 
+    using ResidualNBInfo = typename LocalResidual::ResidualNBInfo;
     struct NeighborInfo
     {
         unsigned int neighbor;
-        double trans;
+        ResidualNBInfo res_nbinfo;
-        double faceArea;
-        FaceDir::DirEnum faceDirection;
         MatrixBlock* matBlockAddress;
     };
     SparseTable<NeighborInfo> neighborInfo_;