energymodule.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_ENERGY_MODULE_HH
#define EWOMS_ENERGY_MODULE_HH

#include <dune/common/fvector.hh>

#include <opm/material/common/Valgrind.hpp>

#include <opm/models/common/multiphasebaseproperties.hh>
#include <opm/models/common/quantitycallbacks.hh>

#include <opm/models/discretization/common/fvbaseproperties.hh>

#include <cassert>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <string>
#include <type_traits>

namespace Opm {
template <class TypeTag, bool enableEnergy>
class EnergyModule;

template <class TypeTag>
class EnergyModule<TypeTag, /*enableEnergy=*/false>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using Model = GetPropType<TypeTag, Properties::Model>;

    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };

public:
    static void registerParameters()
    {}

    static std::string primaryVarName(unsigned)
    { return ""; }

    static std::string eqName(unsigned)
    { return ""; }

    static Scalar primaryVarWeight(const Model&,
                                   unsigned,
                                   unsigned)
    { return -1; }

    static Scalar eqWeight(const Model&,
                           unsigned,
                           unsigned)
    { return -1; }

    template <class FluidState>
    static void setPriVarTemperatures(PrimaryVariables&,
                                      const FluidState&)
    {}

    template <class RateVector, class FluidState>
    static void setEnthalpyRate(RateVector&,
                                const FluidState&,
                                unsigned,
                                const Evaluation&)
    {}

    static void setEnthalpyRate(RateVector&,
                                const Evaluation&)
    {}

    static void addToEnthalpyRate(RateVector&,
                                  const Evaluation&)
    {}

    static Scalar thermalConductionRate(const ExtensiveQuantities&)
    { return 0.0; }

    template <class LhsEval>
    static void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>&,
                                const IntensiveQuantities&,
                                unsigned)
    {}

    template <class LhsEval, class Scv>
    static void addFracturePhaseStorage(Dune::FieldVector<LhsEval, numEq>&,
                                        const IntensiveQuantities&,
                                        const Scv&,
                                        unsigned)
    {}

    template <class LhsEval>
    static void addSolidEnergyStorage(Dune::FieldVector<LhsEval, numEq>&,
                                    const IntensiveQuantities&)
    {}

    template <class Context>
    static void addAdvectiveFlux(RateVector&,
                                 const Context&,
                                 unsigned,
                                 unsigned)
    {}

    template <class Context>
    static void handleFractureFlux(RateVector&,
                                   const Context&,
                                   unsigned,
                                   unsigned)
    {}

    template <class Context>
    static void addDiffusiveFlux(RateVector&,
                                 const Context&,
                                 unsigned,
                                 unsigned)
    {}
};

template <class TypeTag>
class EnergyModule<TypeTag, /*enableEnergy=*/true>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using Model = GetPropType<TypeTag, Properties::Model>;

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

    using Toolbox = Opm::MathToolbox<Evaluation>;

public:
    static void registerParameters()
    {}

    static std::string primaryVarName(unsigned pvIdx)
    {
        if (pvIdx == temperatureIdx) {
            return "temperature";
        }
        return "";
    }

    static std::string eqName(unsigned eqIdx)
    {
        if (eqIdx == energyEqIdx) {
            return "energy";
        }
        return "";
    }

    static Scalar primaryVarWeight(const Model& model, unsigned globalDofIdx, unsigned pvIdx)
    {
        if (pvIdx != temperatureIdx) {
            return -1;
        }

        // make the weight of the temperature primary variable inversly proportional to its value
        return std::max(static_cast<Scalar>(1.0) / 1000,
                        1.0 / model.solution(/*timeIdx=*/0)[globalDofIdx][temperatureIdx]);
    }

    static Scalar eqWeight(const Model&,
                           unsigned,
                           unsigned eqIdx)
    {
        if (eqIdx != energyEqIdx) {
            return -1;
        }

        // approximate change of internal energy of 1kg of liquid water for a temperature
        // change of 30K
        return static_cast<Scalar>(1.0) / (4.184e3 * 30.0);
    }

    static void setEnthalpyRate(RateVector& rateVec, const Evaluation& rate)
    { rateVec[energyEqIdx] = rate; }

    static void addToEnthalpyRate(RateVector& rateVec, const Evaluation& rate)
    { rateVec[energyEqIdx] += rate; }

    static Evaluation thermalConductionRate(const ExtensiveQuantities& extQuants)
    { return -extQuants.temperatureGradNormal() * extQuants.thermalConductivity(); }

    template <class RateVector, class FluidState>
    static void setEnthalpyRate(RateVector& rateVec,
                                const FluidState& fluidState,
                                unsigned phaseIdx,
                                const Evaluation& volume)
    {
        rateVec[energyEqIdx] =
            volume *
            fluidState.density(phaseIdx) *
            fluidState.enthalpy(phaseIdx);
    }

    template <class FluidState>
    static void setPriVarTemperatures(PrimaryVariables& priVars,
                                      const FluidState& fs)
    {
        priVars[temperatureIdx] = Toolbox::value(fs.temperature(/*phaseIdx=*/0));
#ifndef NDEBUG
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            assert(std::abs(Toolbox::value(fs.temperature(/*phaseIdx=*/0))
                            - Toolbox::value(fs.temperature(phaseIdx))) < 1e-30);
        }
#endif
    }

    template <class LhsEval>
    static void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
                                const IntensiveQuantities& intQuants,
                                unsigned phaseIdx)
    {
        const auto& fs = intQuants.fluidState();
        storage[energyEqIdx] +=
            Toolbox::template decay<LhsEval>(fs.density(phaseIdx)) *
            Toolbox::template decay<LhsEval>(fs.internalEnergy(phaseIdx)) *
            Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx)) *
            Toolbox::template decay<LhsEval>(intQuants.porosity());
    }

    template <class Scv, class LhsEval>
    static void addFracturePhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
                                        const IntensiveQuantities& intQuants,
                                        const Scv& scv,
                                        unsigned phaseIdx)
    {
        const auto& fs = intQuants.fractureFluidState();
        storage[energyEqIdx] +=
            Toolbox::template decay<LhsEval>(fs.density(phaseIdx)) *
            Toolbox::template decay<LhsEval>(fs.internalEnergy(phaseIdx)) *
            Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx)) *
            Toolbox::template decay<LhsEval>(intQuants.fracturePorosity()) *
            Toolbox::template decay<LhsEval>(intQuants.fractureVolume()) / scv.volume();
    }

    template <class LhsEval>
    static void addSolidEnergyStorage(Dune::FieldVector<LhsEval, numEq>& storage,
                                    const IntensiveQuantities& intQuants)
    { storage[energyEqIdx] += Opm::decay<LhsEval>(intQuants.solidInternalEnergy()); }

    template <class Context>
    static void addAdvectiveFlux(RateVector& flux,
                                 const Context& context,
                                 unsigned spaceIdx,
                                 unsigned timeIdx)
    {
        const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);

        // advective energy flux in all phases
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!context.model().phaseIsConsidered(phaseIdx)) {
                continue;
            }

            // intensive quantities of the upstream and the downstream DOFs
            const unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
            const IntensiveQuantities& up = context.intensiveQuantities(upIdx, timeIdx);

            flux[energyEqIdx] +=
                extQuants.volumeFlux(phaseIdx) *
                up.fluidState().enthalpy(phaseIdx) *
                up.fluidState().density(phaseIdx);
        }
    }

    template <class Context>
    static void handleFractureFlux(RateVector& flux,
                                   const Context& context,
                                   unsigned spaceIdx,
                                   unsigned timeIdx)
    {
        const auto& scvf = context.stencil(timeIdx).interiorFace(spaceIdx);
        const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);

        // reduce the energy flux in the matrix by the half the width occupied by the
        // fracture
        flux[energyEqIdx] *= 1 - extQuants.fractureWidth() / (2 * scvf.area());

        // advective energy flux in all phases
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!context.model().phaseIsConsidered(phaseIdx)) {
                continue;
            }

            // intensive quantities of the upstream and the downstream DOFs
            const unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
            const IntensiveQuantities& up = context.intensiveQuantities(upIdx, timeIdx);

            flux[energyEqIdx] +=
                extQuants.fractureVolumeFlux(phaseIdx) *
                up.fluidState().enthalpy(phaseIdx) *
                up.fluidState().density(phaseIdx);
        }
    }

    template <class Context>
    static void addDiffusiveFlux(RateVector& flux,
                                 const Context& context,
                                 unsigned spaceIdx,
                                 unsigned timeIdx)
    {
        const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);

        // diffusive energy flux
        flux[energyEqIdx] += -extQuants.temperatureGradNormal() * extQuants.thermalConductivity();
    }
};

template <unsigned PVOffset, bool enableEnergy>
struct EnergyIndices;

template <unsigned PVOffset>
struct EnergyIndices<PVOffset, /*enableEnergy=*/false>
{
    static constexpr unsigned temperatureIdx = std::numeric_limits<unsigned>::max();

    enum { energyEqIdx = -1000 };

protected:
    enum { numEq_ = 0 };
};

template <unsigned PVOffset>
struct EnergyIndices<PVOffset, /*enableEnergy=*/true>
{
    static constexpr unsigned temperatureIdx = PVOffset;

    enum { energyEqIdx = PVOffset };

protected:
    enum { numEq_ = 1 };
};

template <class TypeTag, bool enableEnergy>
class EnergyIntensiveQuantities;

template <class TypeTag>
class EnergyIntensiveQuantities<TypeTag, /*enableEnergy=*/false>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;

    using Toolbox = Opm::MathToolbox<Evaluation>;

public:
    Evaluation solidInternalEnergy() const
    {
        throw std::logic_error("solidInternalEnergy() does not make sense for isothermal models");
    }

    Evaluation thermalConductivity() const
    {
        throw std::logic_error("thermalConductivity() does not make sense for isothermal models");
    }

protected:
    template <class FluidState, class Context>
    static void updateTemperatures_(FluidState& fluidState,
                                    const Context& context,
                                    unsigned spaceIdx,
                                    unsigned timeIdx)
    {
        const Scalar T = context.problem().temperature(context, spaceIdx, timeIdx);
        fluidState.setTemperature(Toolbox::createConstant(T));
    }

    template <class FluidState>
    void update_(FluidState&,
                 typename FluidSystem::template ParameterCache<typename FluidState::ValueType>&,
                 const ElementContext&,
                 unsigned,
                 unsigned)
    {}
};

template <class TypeTag>
class EnergyIntensiveQuantities<TypeTag, /*enableEnergy=*/true>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
    using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;

    enum { numPhases = FluidSystem::numPhases };
    static constexpr unsigned temperatureIdx = Indices::temperatureIdx;

    using Toolbox = Opm::MathToolbox<Evaluation>;

protected:
    template <class FluidState, class Context>
    static void updateTemperatures_(FluidState& fluidState,
                                    const Context& context,
                                    unsigned spaceIdx,
                                    unsigned timeIdx)
    {
        const auto& priVars = context.primaryVars(spaceIdx, timeIdx);
        if constexpr (std::is_same_v<Evaluation, Scalar>) { // finite differences
            fluidState.setTemperature(Toolbox::createConstant(priVars[temperatureIdx]));
        }
        else {
            // automatic differentiation
            if (timeIdx == 0) {
                fluidState.setTemperature(Toolbox::createVariable(priVars[temperatureIdx], temperatureIdx));
            }
            else {
                fluidState.setTemperature(Toolbox::createConstant(priVars[temperatureIdx]));
            }
        }
    }

    template <class FluidState>
    void update_(FluidState& fs,
                 typename FluidSystem::template ParameterCache<typename FluidState::ValueType>& paramCache,
                 const ElementContext& elemCtx,
                 unsigned dofIdx,
                 unsigned timeIdx)
    {
        // set the specific enthalpies of the fluids
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
                continue;
            }

            fs.setEnthalpy(phaseIdx,
                           FluidSystem::enthalpy(fs, paramCache, phaseIdx));
        }

        // compute and set the volumetric internal energy of the solid phase
        const auto& problem = elemCtx.problem();
        const auto& solidEnergyParams = problem.solidEnergyLawParams(elemCtx, dofIdx, timeIdx);
        const auto& thermalCondParams = problem.thermalConductionLawParams(elemCtx, dofIdx, timeIdx);

        solidInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyParams, fs);
        thermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalCondParams, fs);

        Opm::Valgrind::CheckDefined(solidInternalEnergy_);
        Opm::Valgrind::CheckDefined(thermalConductivity_);
    }

public:
    const Evaluation& solidInternalEnergy() const
    { return solidInternalEnergy_; }

    const Evaluation& thermalConductivity() const
    { return thermalConductivity_; }

private:
    Evaluation solidInternalEnergy_;
    Evaluation thermalConductivity_;
};

template <class TypeTag, bool enableEnergy>
class EnergyExtensiveQuantities;

template <class TypeTag>
class EnergyExtensiveQuantities<TypeTag, /*enableEnergy=*/false>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;

protected:
    void update_(const ElementContext&,
                 unsigned,
                 unsigned)
    {}

    template <class Context, class FluidState>
    void updateBoundary_(const Context&,
                         unsigned,
                         unsigned,
                         const FluidState&)
    {}

public:
    Scalar temperatureGradNormal() const
    {
        throw std::logic_error("Calling temperatureGradNormal() does not make sense "
                               "for isothermal models");
    }

    Scalar thermalConductivity() const
    {
        throw std::logic_error("Calling thermalConductivity() does not make sense for "
                               "isothermal models");
    }
};

template <class TypeTag>
class EnergyExtensiveQuantities<TypeTag, /*enableEnergy=*/true>
{
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;

    enum { dimWorld = GridView::dimensionworld };
    using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
    using DimVector = Dune::FieldVector<Scalar, dimWorld>;

protected:
    void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
    {
        const auto& gradCalc = elemCtx.gradientCalculator();
        TemperatureCallback<TypeTag> temperatureCallback(elemCtx);

        EvalDimVector temperatureGrad;
        gradCalc.calculateGradient(temperatureGrad,
                                   elemCtx,
                                   faceIdx,
                                   temperatureCallback);

        // scalar product of temperature gradient and scvf normal
        const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(faceIdx);

        temperatureGradNormal_ = 0;
        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            temperatureGradNormal_ += (face.normal()[dimIdx]*temperatureGrad[dimIdx]);
        }

        const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
        const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
        const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);

        // arithmetic mean
        thermalConductivity_ = 0.5 * (intQuantsInside.thermalConductivity() +
                                      intQuantsOutside.thermalConductivity());
        Opm::Valgrind::CheckDefined(thermalConductivity_);
    }

    template <class Context, class FluidState>
    void updateBoundary_(const Context& context, unsigned bfIdx, unsigned timeIdx, const FluidState& fs)
    {
        const auto& stencil = context.stencil(timeIdx);
        const auto& face = stencil.boundaryFace(bfIdx);

        const auto& elemCtx = context.elementContext();
        const unsigned insideScvIdx = face.interiorIndex();
        const auto& insideScv = elemCtx.stencil(timeIdx).subControlVolume(insideScvIdx);

        const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
        const auto& fsInside = intQuantsInside.fluidState();

        // distance between the center of the SCV and center of the boundary face
        DimVector distVec = face.integrationPos();
        distVec -= insideScv.geometry().center();

        Scalar dist = 0;
        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            dist += distVec[dimIdx] * face.normal()[dimIdx];
        }

        // if the following assertation triggers, the center of the
        // center of the interior SCV was not inside the element!
        assert(dist > 0);

        // calculate the temperature gradient using two-point gradient
        // appoximation
        temperatureGradNormal_ =
            (fs.temperature(/*phaseIdx=*/0) - fsInside.temperature(/*phaseIdx=*/0)) / dist;

        // take the value for thermal conductivity from the interior finite volume
        thermalConductivity_ = intQuantsInside.thermalConductivity();
    }

public:
    const Evaluation& temperatureGradNormal() const
    { return temperatureGradNormal_; }

    const Evaluation& thermalConductivity() const
    { return thermalConductivity_; }

private:
    Evaluation temperatureGradNormal_;
    Evaluation thermalConductivity_;
};

} // namespace Opm

#endif