blackoilenergymodules.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 2 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_BLACK_OIL_ENERGY_MODULE_HH
#define EWOMS_BLACK_OIL_ENERGY_MODULE_HH

#include <dune/common/fvector.hh>

#include <opm/common/ErrorMacros.hpp>
#include <opm/common/utility/gpuDecorators.hpp>

#include <opm/material/common/Tabulated1DFunction.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>

#include <opm/models/blackoil/blackoilmodules.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>
#include <opm/models/common/quantitycallbacks.hh>
#include <opm/models/discretization/common/linearizationtype.hh>
#include <opm/models/io/vtkblackoilenergymodule.hpp>
#include <opm/material/thermal/EnergyModuleType.hpp>

#include <cassert>
#include <cmath>
#include <istream>
#include <memory>
#include <ostream>
#include <stdexcept>
#include <string>


namespace Opm {

template <class TypeTag>
class BlackOilEnergyModule<TypeTag, EnergyModules::FullyImplicitThermal>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using Model = GetPropType<TypeTag, Properties::Model>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;

    static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
    static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;

    static constexpr unsigned enableFullyImplicitThermal = true;
    static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
    static constexpr unsigned numPhases = FluidSystem::numPhases;

public:
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;

    static void registerParameters()
    {
        VtkBlackOilEnergyModule<TypeTag>::registerParameters();
    }

    static void registerOutputModules(Model& model,
                                      Simulator& simulator)
    {
        model.addOutputModule(std::make_unique<VtkBlackOilEnergyModule<TypeTag>>(simulator));
    }

    static OPM_HOST_DEVICE bool primaryVarApplies(unsigned pvIdx)
    {
        return pvIdx == temperatureIdx;
    }

    static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
    {
        assert(primaryVarApplies(pvIdx));

        return "temperature";
    }

    static Scalar primaryVarWeight([[maybe_unused]] unsigned pvIdx)
    {
        assert(primaryVarApplies(pvIdx));

        // TODO: it may be beneficial to chose this differently.
        return static_cast<Scalar>(1.0);
    }

    static OPM_HOST_DEVICE bool eqApplies(unsigned eqIdx)
    {
        return eqIdx == contiEnergyEqIdx;
    }

    static std::string eqName([[maybe_unused]] unsigned eqIdx)
    {
        assert(eqApplies(eqIdx));

        return "conti^energy";
    }

    static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
    {
        assert(eqApplies(eqIdx));

        return 1.0;
    }

    // must be called after water storage is computed
    template <class StorageType>
    OPM_HOST_DEVICE static void addStorage(StorageType& storage,
                                           const IntensiveQuantities& intQuants)
    {
        using LhsEval = typename StorageType::value_type;
        const FluidSystem& fsys = intQuants.getFluidSystem();

        const auto& poro = decay<LhsEval>(intQuants.porosity());

        // accumulate the internal energy of the fluids
        const auto& fs = intQuants.fluidState();
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
            if (!fsys.phaseIsActive(phaseIdx)) {
                continue;
            }

            const auto& u = decay<LhsEval>(fs.internalEnergy(phaseIdx));
            const auto& S = decay<LhsEval>(fs.saturation(phaseIdx));
            const auto& rho = decay<LhsEval>(fs.density(phaseIdx));

            storage[contiEnergyEqIdx] += poro*S*u*rho;
        }

        // add the internal energy of the rock
        const Scalar rockFraction = intQuants.rockFraction();
        const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
        storage[contiEnergyEqIdx] += rockFraction * uRock;
        storage[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
    }

    OPM_HOST_DEVICE static void computeFlux(RateVector& flux,
                                            const ElementContext& elemCtx,
                                            unsigned scvfIdx,
                                            unsigned timeIdx)
    {
        flux[contiEnergyEqIdx] = 0.0;

        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
        const unsigned focusIdx = elemCtx.focusDofIndex();
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!FluidSystem::phaseIsActive(phaseIdx)) {
                continue;
            }

            const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
            if (upIdx == focusIdx) {
                addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
            }
            else {
                addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
            }
        }

        // diffusive energy flux
        flux[contiEnergyEqIdx] += extQuants.energyFlux();
        flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
    }

    template<class RateVectorT>
    OPM_HOST_DEVICE static void addHeatFlux(RateVectorT& flux,
                                            const Evaluation& heatFlux)
    {
        // diffusive energy flux
        flux[contiEnergyEqIdx] += heatFlux;
        flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
    }

    template <class UpEval, class RateVectorT, class Eval, class FluidState>
    OPM_HOST_DEVICE static void addPhaseEnthalpyFluxes_(RateVectorT& 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>
    OPM_HOST_DEVICE static void addPhaseEnthalpyFlux_(RateVector& flux,
                                                      unsigned phaseIdx,
                                                      const ElementContext& elemCtx,
                                                      unsigned scvfIdx,
                                                      unsigned timeIdx)
    {
        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
        const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
        const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
        const auto& fs = up.fluidState();
        const auto& volFlux = extQuants.volumeFlux(phaseIdx);
        addPhaseEnthalpyFluxes_<UpstreamEval>(flux,
                                              phaseIdx,
                                              volFlux,
                                              fs);
    }

    OPM_HOST_DEVICE static void addToEnthalpyRate(RateVector& flux,
                                                  const Evaluation& hRate)
    {
        flux[contiEnergyEqIdx] += hRate;
    }

    template <class FluidState>
    OPM_HOST_DEVICE static void assignPrimaryVars(PrimaryVariables& priVars,
                                                  const FluidState& fluidState)
    {
        priVars[temperatureIdx] = getValue(fluidState.temperature(/*phaseIdx=*/0));
    }

    OPM_HOST_DEVICE static void updatePrimaryVars(PrimaryVariables& newPv,
                                                  const PrimaryVariables& oldPv,
                                                  const EqVector& delta)
    {
        // do a plain unchopped Newton update
        newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
    }

    OPM_HOST_DEVICE static Scalar computeUpdateError(const PrimaryVariables&,
                                                     const EqVector&)
    {
        // do not consider consider the cange of energy primary variables for
        // convergence
        // TODO: maybe this should be changed
        return static_cast<Scalar>(0.0);
    }

    OPM_HOST_DEVICE static Scalar computeResidualError(const EqVector& resid)
    {
        // do not weight the residual of energy when it comes to convergence
        return std::abs(scalarValue(resid[contiEnergyEqIdx]));
    }

    template <class DofEntity>
    static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
    {
        const unsigned dofIdx = model.dofMapper().index(dof);
        const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
        outstream << priVars[temperatureIdx];
    }

    template <class DofEntity>
    static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
    {
        const unsigned dofIdx = model.dofMapper().index(dof);
        PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
        PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];

        instream >> priVars0[temperatureIdx];

        // set the primary variables for the beginning of the current time step.
        priVars1 = priVars0[temperatureIdx];
    }
};

template <class TypeTag>
class BlackOilEnergyIntensiveQuantities<TypeTag, EnergyModules::FullyImplicitThermal>
{
    using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;

    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
    using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;

    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
    static constexpr unsigned temperatureIdx = Indices::temperatureIdx;

public:
    BlackOilEnergyIntensiveQuantities(Evaluation rockInternalEnergy,
                                      Evaluation totalThermalConductivity,
                                      Scalar rockFraction)
        : rockInternalEnergy_(rockInternalEnergy)
        , totalThermalConductivity_(totalThermalConductivity)
        , rockFraction_(rockFraction)
    {
    }

    template <class OtherTypeTag>
    BlackOilEnergyIntensiveQuantities(
        const BlackOilEnergyIntensiveQuantities<OtherTypeTag, EnergyModules::FullyImplicitThermal>& other)
        : rockInternalEnergy_(other.rockInternalEnergy())
        , totalThermalConductivity_(other.totalThermalConductivity())
        , rockFraction_(other.rockFraction())
    {}

    BlackOilEnergyIntensiveQuantities() = default;

    OPM_HOST_DEVICE void updateTemperature_(const ElementContext& elemCtx,
                                            unsigned dofIdx,
                                            unsigned timeIdx)
    {
        auto& fs = asImp_().fluidState_;
        const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);

        // set temperature
        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, elemCtx.linearizationType()));
    }

    OPM_HOST_DEVICE 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));
    }

    OPM_HOST_DEVICE void updateEnergyQuantities_(const ElementContext& elemCtx,
                                                 unsigned dofIdx,
                                                 unsigned timeIdx)
    {
        updateEnergyQuantities_(elemCtx.problem(), elemCtx.globalSpaceIndex(dofIdx, timeIdx), timeIdx);
    }

    OPM_HOST_DEVICE 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 (!asImp_().getFluidSystem().phaseIsActive(phaseIdx)) {
                continue;
            }

            const auto& h = asImp_().getFluidSystem().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);
    }

    OPM_HOST_DEVICE const Evaluation& rockInternalEnergy() const
    { return rockInternalEnergy_; }

    OPM_HOST_DEVICE const Evaluation& totalThermalConductivity() const
    { return totalThermalConductivity_; }

    OPM_HOST_DEVICE Scalar rockFraction() const
    { return rockFraction_; }

protected:
    OPM_HOST_DEVICE Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }

    Evaluation rockInternalEnergy_;
    Evaluation totalThermalConductivity_;
    Scalar rockFraction_;
};

template <class TypeTag>
class BlackOilEnergyIntensiveQuantities<TypeTag, EnergyModules::ConstantTemperature>
{
    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>;

public:

    OPM_HOST_DEVICE void updateTemperature_(const ElementContext& elemCtx,
                                            unsigned dofIdx,
                                            unsigned timeIdx)
    {
        updateTemperature_(elemCtx.problem(), elemCtx.globalSpaceIndex(dofIdx, timeIdx), timeIdx);
    }

    template<class Problem>
    OPM_HOST_DEVICE void updateTemperature_(const Problem& problem,
                                            [[maybe_unused]] const PrimaryVariables& priVars,
                                            unsigned globalDofIdx,
                                            unsigned timeIdx,
                                            [[maybe_unused]] const LinearizationType& lintype)
    {
        updateTemperature_(problem, globalDofIdx, timeIdx);
    }

    OPM_HOST_DEVICE void updateTemperature_(const Problem& problem,
                                            unsigned globalDofIdx,
                                            unsigned timeIdx)
    {
        auto& fs = asImp_().fluidState_;
        const Scalar T = problem.temperature(globalDofIdx, timeIdx);
        fs.setTemperature(T);
    }

    OPM_HOST_DEVICE void updateEnergyQuantities_(const ElementContext&,
                                                 unsigned,
                                                 unsigned,
                                                 const typename FluidSystem::template ParameterCache<Evaluation>&)
    {}

    OPM_HOST_DEVICE const Evaluation& rockInternalEnergy() const
    {
        OPM_THROW(std::logic_error,
                  "Requested the rock internal energy, which is "
                  "unavailable because energy is not conserved");
    }

    OPM_HOST_DEVICE const Evaluation& totalThermalConductivity() const
    {
        OPM_THROW(std::logic_error,
                  "Requested the total thermal conductivity, which is "
                  "unavailable because energy is not conserved");
    }

protected:
    OPM_HOST_DEVICE Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }
};

template <class TypeTag>
class BlackOilEnergyIntensiveQuantities<TypeTag, EnergyModules::SequentialImplicitThermal>
{
    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 SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
    using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;

public:

    OPM_HOST_DEVICE void updateTemperature_(const Problem& problem,
                                            unsigned globalDofIdx,
                                            unsigned timeIdx)
    {
        // update the temperature for output (without derivatives)
        auto& fs = asImp_().fluidState_;
        fs.setTemperature(problem.temperature(globalDofIdx, timeIdx));
    }

    OPM_HOST_DEVICE void updateTemperature_(const ElementContext& elemCtx,
                                            unsigned dofIdx,
                                            unsigned timeIdx)
    {
        updateTemperature_(elemCtx.problem(), elemCtx.globalSpaceIndex(dofIdx, timeIdx), timeIdx);
    }

    template<class Problem>
    OPM_HOST_DEVICE void updateTemperature_(const Problem& problem,
                                            [[maybe_unused]] const PrimaryVariables& priVars,
                                            unsigned globalDofIdx,
                                            unsigned timeIdx,
                                            [[maybe_unused]] const LinearizationType& lintype)
    {
        updateTemperature_(problem, globalDofIdx, timeIdx);
    }

    OPM_HOST_DEVICE void updateEnergyQuantities_([[maybe_unused]] const ElementContext& elemCtx,
                                                 [[maybe_unused]] unsigned dofIdx,
                                                 [[maybe_unused]] unsigned timeIdx)
    {
    }

    OPM_HOST_DEVICE void updateEnergyQuantities_([[maybe_unused]] const Problem& problem,
                                                 [[maybe_unused]] const unsigned globalSpaceIdx,
                                                 [[maybe_unused]] const unsigned timeIdx)
    {
    }

    OPM_HOST_DEVICE const Evaluation& rockInternalEnergy() const
    {
        OPM_THROW(std::logic_error,
                  "Requested the rock internal energy, which is "
                  "unavailable because energy is not conserved");
    }

    OPM_HOST_DEVICE const Evaluation& totalThermalConductivity() const
    {
        OPM_THROW(std::logic_error,
                  "Requested the total thermal conductivity, which is "
                  "unavailable because energy is not conserved");
    }

protected:
    OPM_HOST_DEVICE Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }

};

template <class TypeTag>
class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::FullyImplicitThermal>
{
    using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;

    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;

public:
    template<class Evaluation, class FluidState, class IntensiveQuantities>
    OPM_HOST_DEVICE 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)
    {
        Evaluation deltaT;
        if (focusDofIndex == inIdx) {
            deltaT = decay<Scalar>(exFs.temperature(/*phaseIdx=*/0)) -
                     inFs.temperature(/*phaseIdx=*/0);
        }
        else if (focusDofIndex == exIdx) {
            deltaT = exFs.temperature(/*phaseIdx=*/0) -
                     decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
        }
        else {
            deltaT = decay<Scalar>(exFs.temperature(/*phaseIdx=*/0)) -
                     decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
        }

        Evaluation inLambda;
        if (focusDofIndex == inIdx) {
            inLambda = inIq.totalThermalConductivity();
        }
        else {
            inLambda = decay<Scalar>(inIq.totalThermalConductivity());
        }

        Evaluation exLambda;
        if (focusDofIndex == exIdx) {
            exLambda = exIq.totalThermalConductivity();
        }
        else {
            exLambda = decay<Scalar>(exIq.totalThermalConductivity());
        }

        Evaluation H;
        const Evaluation& inH = inLambda*inAlpha;
        const Evaluation& exH = exLambda*outAlpha;
        if (inH > 0 && exH > 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.
            H = 1.0 / (1.0 / inH + 1.0 / exH);
        }
        else {
            H = 0.0;
        }

        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>
    void updateEnergyBoundary(const Context& ctx,
                              unsigned scvfIdx,
                              unsigned timeIdx,
                              const BoundaryFluidState& boundaryFs)
    {
        const auto& stencil = ctx.stencil(timeIdx);
        const auto& scvf = stencil.boundaryFace(scvfIdx);

        const unsigned inIdx = scvf.interiorIndex();
        const auto& inIq = ctx.intensiveQuantities(inIdx, timeIdx);
        const auto& focusDofIdx = ctx.focusDofIndex();
        const Scalar alpha = ctx.problem().thermalHalfTransmissibilityBoundary(ctx, scvfIdx);
        updateEnergyBoundary(energyFlux_, inIq, focusDofIdx, inIdx, alpha, boundaryFs);
    }

    template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
    OPM_HOST_DEVICE 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 (focusDofIndex == inIdx) {
            deltaT = boundaryFs.temperature(/*phaseIdx=*/0) -
                     inFs.temperature(/*phaseIdx=*/0);
        }
        else {
            deltaT = decay<Scalar>(boundaryFs.temperature(/*phaseIdx=*/0)) -
                     decay<Scalar>(inFs.temperature(/*phaseIdx=*/0));
        }

        Evaluation lambda;
        if (focusDofIndex == inIdx) {
            lambda = inIq.totalThermalConductivity();
        }
        else {
            lambda = decay<Scalar>(inIq.totalThermalConductivity());
        }

        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.
            energyFlux = deltaT * lambda * -alpha;
        }
        else {
            energyFlux = 0.0;
        }
    }

    const Evaluation& energyFlux()  const
    { return energyFlux_; }

private:
    Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }

    Evaluation energyFlux_;
};

template <class TypeTag>
class BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::SequentialImplicitThermal>
{
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;

public:
    template<class Evaluation, class FluidState, class IntensiveQuantities>
    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)
    {
        // Uses the same caculations as the FullyImplicitThermal approach
        BlackOilEnergyExtensiveQuantities<TypeTag, EnergyModules::FullyImplicitThermal>::updateEnergy(energyFlux,
                     focusDofIndex,
                     inIdx,
                     exIdx,
                     inIq,
                     exIq,
                     inFs,
                     exFs,
                     inAlpha,
                     outAlpha,
                     faceArea);
    }

    void updateEnergy(const ElementContext&,
                      unsigned,
                      unsigned)
    { } // Old interface still used output code for fluxes. But energy flux is not used. i.e. do nothing

    template <class Context, class BoundaryFluidState>
    void updateEnergyBoundary(const Context&,
                              unsigned,
                              unsigned,
                              const BoundaryFluidState&)
    { }

    template <class Evaluation, class BoundaryFluidState, class IntensiveQuantities>
    static void updateEnergyBoundary(Evaluation& /*heatFlux*/,
                                     const IntensiveQuantities& /*inIq*/,
                                     unsigned /*focusDofIndex*/,
                                     unsigned /*inIdx*/,
                                     unsigned /*timeIdx*/,
                                     Scalar /*alpha*/,
                                     const BoundaryFluidState& /*boundaryFs*/)
    { }

    OPM_HOST_DEVICE const Evaluation& energyFlux() const
    { OPM_THROW(std::logic_error, "Requested the energy flux, but energy is not conserved"); }
};

} // namespace Opm

#endif