implement the ebos part of energy conservation

I.e. everything which is ECL specific.

ebos/eclfluxmodule.hh

@@ -115,6 +115,7 @@ class EclTransExtensiveQuantities
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
     enum { numPhases = FluidSystem::numPhases };
     enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
 
     typedef Opm::MathToolbox<Evaluation> Toolbox;
     typedef Dune::FieldVector<Scalar, dimWorld> DimVector;

ebos/eclpeacemanwell.hh

@@ -116,8 +116,13 @@ class EclPeacemanWell : public BaseAuxiliaryModule<TypeTag>
 
     static const unsigned numModelEq = GET_PROP_VALUE(TypeTag, NumEq);
     static const unsigned conti0EqIdx = GET_PROP_TYPE(TypeTag, Indices)::conti0EqIdx;
+    static const unsigned contiEnergyEqIdx = GET_PROP_TYPE(TypeTag, Indices)::contiEnergyEqIdx;
 
-    typedef Opm::CompositionalFluidState<Scalar, FluidSystem, /*storeEnthalpy=*/false> FluidState;
+    static constexpr unsigned historySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize);
+
+    static constexpr bool enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy);
+
+    typedef Opm::CompositionalFluidState<Scalar, FluidSystem, /*storeEnthalpy=*/true> FluidState;
     typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
 
     // all quantities that need to be stored per degree of freedom that intersects the
@@ -318,6 +323,17 @@ public:
 
         int wellGlobalDof = AuxModule::localToGlobalDof(/*localDofIdx=*/0);
         sol[wellGlobalDof] = 0.0;
+
+        // make valgrind shut up about the DOFs for the well even if the PrimaryVariables
+        // class contains some "holes" due to alignment
+        Opm::Valgrind::SetDefined(sol[wellGlobalDof]);
+
+        // also apply the initial solution of the well to the "old" time steps
+        for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx) {
+            auto& oldSol = const_cast<SolutionVector&>(simulator_.model().solution(timeIdx));
+
+            oldSol[wellGlobalDof] = sol[wellGlobalDof];
+        }
     }
 
     /*!
@@ -708,6 +724,12 @@ public:
     void setControlMode(ControlMode controlMode)
     { controlMode_ = controlMode; }
 
+    /*!
+     * \brief Set the temperature of the injected fluids [K]
+     */
+    void setTemperature(Scalar value)
+    { wellTemperature_ = value; }
+
     /*!
      * \brief Set the connection transmissibility factor for a given degree of freedom.
      */
@@ -1104,15 +1126,59 @@ public:
         computeVolumetricDofRates_(volumetricRates, actualBottomHolePressure_, tmp);
 
         // convert to mass rates
-        RateVector modelRate;
+        RateVector modelRate(0.0);
         const auto& intQuants = context.intensiveQuantities(dofIdx, timeIdx);
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx))
                 continue;
 
+            // energy is disabled or we have production for the given phase, i.e., we
+            // can use the intensive quantities' fluid state
             modelRate.setVolumetricRate(intQuants.fluidState(), phaseIdx, volumetricRates[phaseIdx]);
-            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
+
-                q[conti0EqIdx + compIdx] += modelRate[conti0EqIdx + compIdx];
+            if (enableEnergy) {
+                if (volumetricRates[phaseIdx] < 0.0) {
+                    // producer
+                    const auto& fs = intQuants.fluidState();
+                    modelRate[contiEnergyEqIdx] += volumetricRates[phaseIdx]*fs.density(phaseIdx)*fs.enthalpy(phaseIdx);
+                }
+                else if (volumetricRates[phaseIdx] > 0.0
+                         && injectedPhaseIdx_ == phaseIdx)
+                {
+                    // injector for the right phase. we need to use the thermodynamic
+                    // quantities from the borehole as upstream
+                    //
+                    // TODO: This is not implemented in a very efficient way, the
+                    // required quantities could be precomputed at initialization!
+                    auto fs = injectionFluidState_;
+
+                    // TODO: maybe we need to use a depth dependent pressure here. the
+                    // difference is probably not very large, and for wells that span
+                    // multiple perforations it is unclear what "well temperature" means
+                    // anyway.
+                    fs.setPressure(phaseIdx, actualBottomHolePressure_);
+
+                    fs.setTemperature(wellTemperature_);
+
+                    typename FluidSystem::template ParameterCache<Evaluation> paramCache;
+                    unsigned globalSpaceIdx = context.globalSpaceIndex(dofIdx, timeIdx);
+                    unsigned pvtRegionIdx = context.primaryVars(dofIdx, timeIdx).pvtRegionIndex();
+                    paramCache.setRegionIndex(pvtRegionIdx);
+                    paramCache.setMaxOilSat(context.problem().maxOilSaturation(globalSpaceIdx));
+                    paramCache.updatePhase(fs, phaseIdx);
+
+                    const auto& rho = FluidSystem::density(fs, paramCache, phaseIdx);
+                    fs.setDensity(phaseIdx, rho);
+
+                    const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
+                    fs.setEnthalpy(phaseIdx, h);
+
+                    modelRate[contiEnergyEqIdx] += volumetricRates[phaseIdx]*fs.density(phaseIdx)*fs.enthalpy(phaseIdx);
+                }
+            }
+
+            for (unsigned eqIdx = 0; eqIdx < modelRate.size(); ++eqIdx)
+                q[conti0EqIdx + eqIdx] += modelRate[conti0EqIdx + eqIdx];
         }
 
         Opm::Valgrind::CheckDefined(q);
@@ -1562,6 +1628,9 @@ protected:
     // the sum of the total volumes of all the degrees of freedoms that interact with the well
     Scalar wellTotalVolume_;
 
+    // the temperature assumed for the fluid (in the case of an injector well)
+    Scalar wellTemperature_;
+
     // The assumed bottom hole and tubing head pressures as specified by the user
     Scalar bhpLimit_;
     Scalar thpLimit_;

ebos/eclproblem.hh

@@ -65,6 +65,8 @@
 #include <ewoms/disc/ecfv/ecfvdiscretization.hh>
 
 #include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
+#include <opm/material/thermal/EclThermalLawManager.hpp>
+
 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 #include <opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp>
@@ -80,6 +82,7 @@
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/parser/eclipse/EclipseState/SummaryConfig/SummaryConfig.hpp>
 #include <opm/material/common/Exceptions.hpp>
+#include <opm/material/common/ConditionalStorage.hpp>
 
 #include <dune/common/version.hh>
 #include <dune/common/fvector.hh>
@@ -129,7 +132,7 @@ SET_TAG_PROP(EclBaseProblem, SpatialDiscretizationSplice, EcfvDiscretization);
 //! for ebos, use automatic differentiation to linearize the system of PDEs
 SET_TAG_PROP(EclBaseProblem, LocalLinearizerSplice, AutoDiffLocalLinearizer);
 
-// Set the material Law
+// Set the material law for fluid fluxes
 SET_PROP(EclBaseProblem, MaterialLaw)
 {
 private:
@@ -147,6 +150,32 @@ public:
     typedef typename EclMaterialLawManager::MaterialLaw type;
 };
 
+// Set the material law for heat heat storage in rock
+SET_PROP(EclBaseProblem, SolidEnergyLaw)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+    typedef Opm::EclThermalLawManager<Scalar, FluidSystem> EclThermalLawManager;
+
+    typedef typename EclThermalLawManager::SolidEnergyLaw type;
+};
+
+// Set the material law for heat conduction
+SET_PROP(EclBaseProblem, ThermalConductionLaw)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+    typedef Opm::EclThermalLawManager<Scalar, FluidSystem> EclThermalLawManager;
+
+    typedef typename EclThermalLawManager::ThermalConductionLaw type;
+};
+
 // ebos can use a slightly faster stencil class because it does not need the normals and
 // the integration points of intersections
 SET_PROP(EclBaseProblem, Stencil)
@@ -267,6 +296,7 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
     enum { numComponents = FluidSystem::numComponents };
     enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
     enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
+    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
     enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
     enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
@@ -281,13 +311,15 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
     typedef typename GridView::template Codim<0>::Entity Element;
     typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
     typedef typename GET_PROP(TypeTag, MaterialLaw)::EclMaterialLawManager EclMaterialLawManager;
-    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+    typedef typename GET_PROP(TypeTag, SolidEnergyLaw)::EclThermalLawManager EclThermalLawManager;
+    typedef typename EclMaterialLawManager::MaterialLawParams MaterialLawParams;
+    typedef typename EclThermalLawManager::SolidEnergyLawParams SolidEnergyLawParams;
+    typedef typename EclThermalLawManager::ThermalConductionLawParams ThermalConductionLawParams;
     typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
-    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
     typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
     typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
 
-
     typedef BlackOilSolventModule<TypeTag> SolventModule;
     typedef BlackOilPolymerModule<TypeTag> PolymerModule;
 
@@ -337,13 +369,14 @@ public:
         , transmissibilities_(simulator.vanguard())
         , thresholdPressures_(simulator)
         , wellManager_(simulator)
+        , eclWriter_(EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput)
+                     ? new EclWriterType(simulator) : nullptr)
         , pffDofData_(simulator.gridView(), this->elementMapper())
     {
         // Tell the black-oil extensions to initialize their internal data structures
         const auto& vanguard = simulator.vanguard();
         SolventModule::initFromDeck(vanguard.deck(), vanguard.eclState());
         PolymerModule::initFromDeck(vanguard.deck(), vanguard.eclState());
-
         if (EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput))
             // create the ECL writer
             eclWriter_.reset(new EclWriterType(simulator));
@@ -418,6 +451,7 @@ public:
         updateElementDepths_();
         readRockParameters_();
         readMaterialParameters_();
+        readThermalParameters_();
         transmissibilities_.finishInit();
         readInitialCondition_();
 
@@ -723,7 +757,7 @@ public:
     { return transmissibilities_.permeability(globalElemIdx); }
 
     /*!
-     * \copydoc BlackOilBaseProblem::transmissibility
+     * \copydoc EclTransmissiblity::transmissibility
      */
     template <class Context>
     Scalar transmissibility(const Context& context,
@@ -734,6 +768,30 @@ public:
         return pffDofData_.get(context.element(), toDofLocalIdx).transmissibility;
     }
 
+    /*!
+     * \copydoc EclTransmissiblity::thermalHalfTransmissibility
+     */
+    template <class Context>
+    Scalar thermalHalfTransmissibility(const Context& context,
+                                       unsigned OPM_OPTIM_UNUSED fromDofLocalIdx,
+                                       unsigned toDofLocalIdx) const
+    {
+        assert(fromDofLocalIdx == 0);
+        return *pffDofData_.get(context.element(), toDofLocalIdx).thermalHalfTrans;
+    }
+
+    /*!
+     * \copydoc EclTransmissiblity::thermalHalfTransmissibility
+     */
+    template <class Context>
+    Scalar thermalHalfTransmissibilityBoundary(const Context& context,
+                                               unsigned dofIdx,
+                                               unsigned boundaryFaceIdx) const
+    {
+        unsigned elemIdx = context.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
+        return transmissibilities_.thermalHalfTransBoundary(elemIdx, boundaryFaceIdx);
+    }
+
     /*!
      * \brief Return a reference to the object that handles the "raw" transmissibilities.
      */
@@ -829,6 +887,31 @@ public:
     const MaterialLawParams& materialLawParams(unsigned globalDofIdx) const
     { return materialLawManager_->materialLawParams(globalDofIdx); }
 
+    /*!
+     * \brief Return the parameters for the heat storage law of the rock
+     */
+    template <class Context>
+    const SolidEnergyLawParams&
+    solidEnergyLawParams(const Context& context OPM_UNUSED,
+                         unsigned spaceIdx OPM_UNUSED,
+                         unsigned timeIdx OPM_UNUSED) const
+    {
+        unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
+        return thermalLawManager_->solidEnergyLawParams(globalSpaceIdx);
+    }
+
+
+    /*!
+     * \copydoc FvBaseMultiPhaseProblem::thermalConductionParams
+     */
+    template <class Context>
+    const ThermalConductionLawParams &
+    thermalConductionLawParams(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
+    {
+        unsigned globalSpaceIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
+        return thermalLawManager_->thermalConductionLawParams(globalSpaceIdx);
+    }
+
     /*!
      * \brief Returns the ECL material law manager
      *
@@ -844,8 +927,6 @@ public:
     std::shared_ptr<EclMaterialLawManager> materialLawManager()
     { return materialLawManager_; }
 
-
-
     /*!
      * \brief Returns the initial solvent saturation for a given a cell index
      */
@@ -1018,6 +1099,8 @@ public:
 
         if (enablePolymer)
              values[Indices::polymerConcentrationIdx] = polymerConcentration_[globalDofIdx];
+
+        values.checkDefined();
     }
 
     /*!
@@ -1076,8 +1159,12 @@ public:
             // convert the source term from the total mass rate of the
             // cell to the one per unit of volume as used by the model.
             unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
-            for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
                 rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
+
+                Opm::Valgrind::CheckDefined(rate[eqIdx]);
+                assert(Opm::isfinite(rate[eqIdx]));
+            }
         }
     }
 
@@ -1374,6 +1461,25 @@ private:
         ////////////////////////////////
     }
 
+    void readThermalParameters_()
+    {
+        if (!enableEnergy)
+            return;
+
+        const auto& vanguard = this->simulator().vanguard();
+        const auto& deck = vanguard.deck();
+        const auto& eclState = vanguard.eclState();
+
+        // fluid-matrix interactions (saturation functions; relperm/capillary pressure)
+        size_t numDof = this->model().numGridDof();
+        std::vector<int> compressedToCartesianElemIdx(numDof);
+        for (unsigned elemIdx = 0; elemIdx < numDof; ++elemIdx)
+            compressedToCartesianElemIdx[elemIdx] = vanguard.cartesianIndex(elemIdx);
+
+        thermalLawManager_ = std::make_shared<EclThermalLawManager>();
+        thermalLawManager_->initFromDeck(deck, eclState, compressedToCartesianElemIdx);
+    }
+
     void updatePorosity_()
     {
         const auto& vanguard = this->simulator().vanguard();
@@ -1397,8 +1503,10 @@ private:
             Scalar poreVolume = porvData[cartElemIdx];
 
             // sum up the pore volume of the active cell and all inactive ones above it
-            // which were disabled due to their pore volume being too small
+            // which were disabled due to their pore volume being too small. If energy is
-            if (eclGrid.getMinpvMode() == Opm::MinpvMode::ModeEnum::OpmFIL) {
+            // conserved, cells are not disabled due to a too small pore volume because
+            // such cells still store and conduct energy.
+            if (!enableEnergy && eclGrid.getMinpvMode() == Opm::MinpvMode::ModeEnum::OpmFIL) {
                 Scalar minPvValue = eclGrid.getMinpvValue();
                 for (int aboveElemCartIdx = static_cast<int>(cartElemIdx) - nx*ny;
                      aboveElemCartIdx >= 0;
@@ -1780,6 +1888,7 @@ private:
 
     struct PffDofData_
     {
+        Opm::ConditionalStorage<enableEnergy, Scalar> thermalHalfTrans;
         Scalar transmissibility;
     };
 
@@ -1798,6 +1907,9 @@ private:
             if (localDofIdx != 0) {
                 unsigned globalCenterElemIdx = elementMapper.index(stencil.entity(/*dofIdx=*/0));
                 dofData.transmissibility = transmissibilities_.transmissibility(globalCenterElemIdx, globalElemIdx);
+
+                if (enableEnergy)
+                    *dofData.thermalHalfTrans = transmissibilities_.thermalHalfTrans(globalCenterElemIdx, globalElemIdx);
             }
         };
 
@@ -1809,6 +1921,7 @@ private:
     EclTransmissibility<TypeTag> transmissibilities_;
 
     std::shared_ptr<EclMaterialLawManager> materialLawManager_;
+    std::shared_ptr<EclThermalLawManager> thermalLawManager_;
 
     EclThresholdPressure<TypeTag> thresholdPressures_;
 
@@ -1833,7 +1946,6 @@ private:
     std::vector<Scalar> lastRs_;
     Scalar maxDRsDt_;
     Scalar maxDRs_;
-
     bool drvdtActive_; // if no, VAPPARS *might* be active
     std::vector<Scalar> lastRv_;
     Scalar maxDRvDt_;

ebos/ecltransmissibility.hh

@@ -28,7 +28,6 @@
 #ifndef EWOMS_ECL_TRANSMISSIBILITY_HH
 #define EWOMS_ECL_TRANSMISSIBILITY_HH
 
-
 #include <ewoms/common/propertysystem.hh>
 
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -39,6 +38,7 @@
 #include <opm/grid/CpGrid.hpp>
 
 #include <opm/material/common/Exceptions.hpp>
+#include <opm/material/common/ConditionalStorage.hpp>
 
 #include <dune/grid/common/mcmgmapper.hh>
 
@@ -57,6 +57,7 @@ NEW_PROP_TAG(Vanguard);
 NEW_PROP_TAG(Grid);
 NEW_PROP_TAG(GridView);
 NEW_PROP_TAG(ElementMapper);
+NEW_PROP_TAG(EnableEnergy);
 }
 
 /*!
@@ -75,6 +76,8 @@ class EclTransmissibility
     typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
     typedef typename GridView::Intersection Intersection;
 
+    static const bool enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy);
+
     // Grid and world dimension
     enum { dimWorld = GridView::dimensionworld };
 
@@ -104,6 +107,11 @@ public:
     { update(); }
 
 
+    /*!
+     * \brief Compute all transmissibilities
+     *
+     * Also, this updates the "thermal half transmissibilities" if energy is enabled.
+     */
     void update()
     {
         const auto& gridView = vanguard_.gridView();
@@ -154,10 +162,18 @@ public:
         trans_.clear();
         trans_.reserve(numElements*3*1.05);
 
+        // if energy is enabled, let's do the same for the "thermal half transmissibilities"
+        if (enableEnergy) {
+            thermalHalfTrans_->clear();
+            thermalHalfTrans_->reserve(numElements*3*1.05);
+        }
+
         // compute the transmissibilities for all intersections
         elemIt = gridView.template begin</*codim=*/ 0>();
         for (; elemIt != elemEndIt; ++elemIt) {
             const auto& elem = *elemIt;
+            unsigned elemIdx = elemMapper.index(elem);
+
             auto isIt = gridView.ibegin(elem);
             const auto& isEndIt = gridView.iend(elem);
             for (; isIt != isEndIt; ++ isIt) {
@@ -165,22 +181,31 @@ public:
                 const auto& intersection = *isIt;
 
                 // ignore boundary intersections for now (TODO?)
-
+                if (!intersection.neighbor())
-                // continue if no neighbor is present
-                if ( ! intersection.neighbor() )
                     continue;
 
-                const auto& inside = intersection.inside();
+                const auto& outsideElem = intersection.outside();
-                const auto& outside = intersection.outside();
+                unsigned outsideElemIdx = elemMapper.index(outsideElem);
-                unsigned insideElemIdx = elemMapper.index(inside);
-                unsigned outsideElemIdx = elemMapper.index(outside);
 
                 // we only need to calculate a face's transmissibility
                 // once...
-                if (insideElemIdx > outsideElemIdx)
+                if (elemIdx > outsideElemIdx)
                     continue;
 
-                unsigned insideCartElemIdx = cartMapper.cartesianIndex(insideElemIdx);
+                // update the "thermal half transmissibility" for the intersection
+                if (enableEnergy) {
+                    const auto& n = intersection.centerUnitOuterNormal();
+                    Scalar A = intersection.geometry().volume();
+
+                    const auto& inPos = elem.geometry().center();
+                    const auto& outPos = outsideElem.geometry().center();
+                    const auto& d = outPos - inPos;
+
+                    (*thermalHalfTrans_)[isId_(elemIdx, outsideElemIdx)] =
+                        A * std::abs(n*d)/(d*d);
+                }
+
+                unsigned insideCartElemIdx = cartMapper.cartesianIndex(elemIdx);
                 unsigned outsideCartElemIdx = cartMapper.cartesianIndex(outsideElemIdx);
 
                 // local indices of the faces of the inside and
@@ -193,7 +218,7 @@ public:
                 DimVector faceAreaNormal;
 
                 typename std::is_same< Grid, Dune::CpGrid> :: type isCpGrid;
-                computeFaceProperties( intersection, insideElemIdx, insideFaceIdx, outsideElemIdx, outsideFaceIdx,
+                computeFaceProperties( intersection, elemIdx, insideFaceIdx, outsideElemIdx, outsideFaceIdx,
                                        faceCenterInside, faceCenterOutside, faceAreaNormal,
                                        isCpGrid );
 
@@ -205,9 +230,9 @@ public:
                                   insideFaceIdx,
                                   distanceVector_(faceCenterInside,
                                                   intersection.indexInInside(),
-                                                  insideElemIdx,
+                                                  elemIdx,
                                                   axisCentroids),
-                                  permeability_[insideElemIdx]);
+                                  permeability_[elemIdx]);
                 computeHalfTrans_(halfTrans2,
                                   faceAreaNormal,
                                   outsideFaceIdx,
@@ -261,17 +286,40 @@ public:
                                                        outsideCartElemIdx,
                                                        faceDir);
 
-                trans_[isId_(insideElemIdx, outsideElemIdx)] = trans;
+                trans_[isId_(elemIdx, outsideElemIdx)] = trans;
             }
         }
     }
 
+    /*!
+     * \brief Return the permeability for an element.
+     */
     const DimMatrix& permeability(unsigned elemIdx) const
     { return permeability_[elemIdx]; }
 
+    /*!
+     * \brief Return the transmissibility for the intersection between two elements.
+     */
     Scalar transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
     { return trans_.at(isId_(elemIdx1, elemIdx2)); }
 
+    /*!
+     * \brief Return the thermal "half transmissibility" for the intersection between two
+     *        elements.
+     *
+     * The "half transmissibility" features all sub-expressions of the "thermal
+     * transmissibility" which can be precomputed and is not dependent on the current
+     * solution:
+     *
+     * H_t = A* (n*d)/(d*d);
+     *
+     * where A is the area of the intersection between the inside and outside elements, n
+     * is the outer unit normal and d is the distance between the centers of the adjacent
+     * elements
+     */
+    Scalar thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
+    { return thermalHalfTrans_->at(isId_(insideElemIdx, outsideElemIdx)); }
+
 private:
     template <class Intersection>
     void computeFaceProperties( const Intersection& intersection,
@@ -484,15 +532,13 @@ private:
             }
             averageNtg[cartesianCellIdx] = ntgCellVolume / totalCellVolume;
         }
-
     }
 
-
-
-
     const Vanguard& vanguard_;
     std::vector<DimMatrix> permeability_;
     std::unordered_map<std::uint64_t, Scalar> trans_;
+    Opm::ConditionalStorage<enableEnergy,
+                            std::unordered_map<std::uint64_t, Scalar> > thermalHalfTrans_;
 };
 
 } // namespace Ewoms

ebos/eclwellmanager.hh

@@ -103,6 +103,9 @@ public:
         {
             const Opm::Well* deckWell = deckSchedule.getWells()[deckWellIdx];
             const std::string& wellName = deckWell->name();
+            Scalar wellTemperature = 273.15 + 15.56; // [K]
+            if (deckWell->isInjector(/*timeStep=*/0))
+                wellTemperature = deckWell->getInjectionProperties(/*timeStep=*/0).temperature;
 
             // set the name of the well but not much else. (i.e., if it is not completed,
             // the well primarily serves as a placeholder.) The big rest of the well is
@@ -110,6 +113,7 @@ public:
             auto well = std::make_shared<Well>(simulator_);
             well->setName(wellName);
             well->setWellStatus(Well::Shut);
+            well->setTemperature(wellTemperature);
 
             wells_.push_back(well);
             wellNameToIndex_[well->name()] = wells_.size() - 1;
@@ -144,6 +148,9 @@ public:
 
             auto well = this->well(deckWell->name());
 
+            if (deckWell->isInjector(episodeIdx))
+                well->setTemperature(deckWell->getInjectionProperties(episodeIdx).temperature);
+
             Opm::WellCommon::StatusEnum deckWellStatus = deckWell->getStatus(episodeIdx);
             switch (deckWellStatus) {
             case Opm::WellCommon::AUTO: