Commit 4919aed4 authored by Timo Koch's avatar Timo Koch
Browse files

[energy] Add concept for energy balance

* Implement fluxvar specializations with energy on
* EnergyLocalResidual helps to add energy contributions efficiently
* The base implicit volVars now exist in an isothermal and a non-isothermal version
* Make 1p non-isothermal cc models run

* TODO
	* Fourier's law for box
	* Adjust volume variables / local residuals of other models
parent 996dfeb4
......@@ -25,29 +25,43 @@
#define DUMUX_DISCRETIZATION_VOLUME_VARIABLES_HH
#include <dumux/implicit/properties.hh>
#include <dumux/common/valgrind.hh>
namespace Dumux
{
namespace Properties
{
NEW_PROP_TAG(FluidSystem);
NEW_PROP_TAG(Indices);
}
// forward declaration
template <class TypeTag, bool enableEnergyBalance>
class ImplicitVolumeVariablesImplementation;
/*!
* \ingroup ImplicitVolumeVariables
* \brief Base class for the model specific class which provides
* access to all volume averaged quantities.
* access to all volume averaged quantities. The volume variables base class
* is specialized for isothermal and non-isothermal models.
*/
template <class TypeTag>
class ImplicitVolumeVariables
{
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
using ImplicitVolumeVariables = ImplicitVolumeVariablesImplementation<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergyBalance)>;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
/*!
* \ingroup ImplicitVolumeVariables
* \brief The isothermal base class
*/
template<class TypeTag>
class ImplicitVolumeVariablesImplementation<TypeTag, false>
{
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
public:
......@@ -72,7 +86,6 @@ public:
const Element &element,
const SubControlVolume &scv)
{
Valgrind::CheckDefined(priVars);
priVars_ = priVars;
extrusionFactor_ = problem.boxExtrusionFactor(element, scv);
}
......@@ -89,9 +102,7 @@ public:
* \param pvIdx The index of the primary variable of interest
*/
Scalar priVar(const int pvIdx) const
{
return priVars_[pvIdx];
}
{ return priVars_[pvIdx]; }
/*!
* \brief Return how much the sub-control volume is extruded.
......@@ -105,15 +116,135 @@ public:
Scalar extrusionFactor() const
{ return extrusionFactor_; }
//! The temperature is obtained from the problem as a constant for isothermal models
static Scalar temperature(const PrimaryVariables &priVars,
const Problem& problem,
const Element &element,
const SubControlVolume &scv)
{
return problem.temperatureAtPos(scv.dofPosition());
}
//! The phase enthalpy is zero for isothermal models
//! This is needed for completing the fluid state
template<class FluidState, class ParameterCache>
static Scalar enthalpy(const FluidState& fluidState,
const ParameterCache& paramCache,
const int phaseIdx)
{
return 0;
}
private:
PrimaryVariables priVars_;
Scalar extrusionFactor_;
};
//! The non-isothermal implicit volume variables base class
template <class TypeTag>
class ImplicitVolumeVariablesImplementation<TypeTag, true>
: public ImplicitVolumeVariablesImplementation<TypeTag, false>
{
using ParentType = ImplicitVolumeVariablesImplementation<TypeTag, false>;
using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum { temperatureIdx = Indices::temperatureIdx };
public:
/*!
* \brief Update all quantities for a given control volume
*
* \param priVars A vector containing the primary variables for the control volume
* \param problem The object specifying the problem which ought to
* be simulated
* \param element An element which contains part of the control volume
* \param fvGeometry The finite volume geometry for the element
* \param scvIdx Local index of the sub control volume which is inside the element
*/
void update(const PrimaryVariables &priVars,
const Problem &problem,
const Element &element,
const SubControlVolume &scv)
{
ParentType::update(priVars, problem, element, scv);
solidHeatCapacity_ = problem.spatialParams().solidHeatCapacity(element, scv);
solidDensity_ = problem.spatialParams().solidDensity(element, scv);
solidThermalConductivity_ = problem.spatialParams().solidThermalConductivity(element, scv);
}
/*!
* \brief If running in valgrind this makes sure that all
* quantities in the volume variables are defined.
* \brief Returns the total internal energy of a phase in the
* sub-control volume.
*
* \param phaseIdx The phase index
*/
void checkDefined() const
Scalar internalEnergy(const int phaseIdx) const
{ return asImp_().fluidState().internalEnergy(phaseIdx); }
/*!
* \brief Returns the total enthalpy of a phase in the sub-control
* volume.
*
* \param phaseIdx The phase index
*/
Scalar enthalpy(const int phaseIdx) const
{ return asImp_().fluidState().enthalpy(phaseIdx); }
/*!
* \brief Returns the total heat capacity \f$\mathrm{[J/(kg K)]}\f$ of the rock matrix in
* the sub-control volume.
*/
Scalar solidHeatCapacity() const
{ return solidHeatCapacity_; }
/*!
* \brief Returns the mass density \f$\mathrm{[kg/m^3]}\f$ of the rock matrix in
* the sub-control volume.
*/
Scalar solidDensity() const
{ return solidDensity_; }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of a fluid phase in
* the sub-control volume.
*/
Scalar fluidThermalConductivity(const int phaseIdx) const
{ return FluidSystem::thermalConductivity(asImp_().fluidState(), phaseIdx); }
/*!
* \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the solid phase in
* the sub-control volume.
*/
Scalar solidThermalConductivity() const
{ return solidThermalConductivity_; }
//! The temperature is a primary variable for non-isothermal models
static Scalar temperature(const PrimaryVariables &priVars,
const Problem& problem,
const Element &element,
const SubControlVolume &scv)
{
#if !defined NDEBUG && HAVE_VALGRIND
Valgrind::CheckDefined(priVars_);
#endif
return priVars[temperatureIdx];
}
//! The phase enthalpy is zero for isothermal models
//! This is needed for completing the fluid state
template<class FluidState, class ParameterCache>
static Scalar enthalpy(const FluidState& fluidState,
const ParameterCache& paramCache,
const int phaseIdx)
{
return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
}
protected:
......@@ -122,10 +253,12 @@ protected:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
PrimaryVariables priVars_;
Scalar extrusionFactor_;
private:
Scalar solidHeatCapacity_;
Scalar solidDensity_;
Scalar solidThermalConductivity_;
};
} // end namespace
} // end namespace Dumux
#endif
......@@ -55,6 +55,7 @@ NEW_PROP_TAG(BaseModel); //!< The type of the base class of the model
NEW_PROP_TAG(NumEq); //!< Number of equations in the system of PDEs
NEW_PROP_TAG(BaseLocalResidual); //!< The type of the base class of the local residual
NEW_PROP_TAG(LocalResidual); //!< The type of the local residual function
NEW_PROP_TAG(EnergyLocalResidual); //! The local residual of the energy equation
NEW_PROP_TAG(LocalJacobian); //!< The type of the local jacobian operator
NEW_PROP_TAG(SubControlVolume);//!< The type of the sub control volume
......@@ -92,6 +93,7 @@ NEW_PROP_TAG(AdvectionType); //! The type for the calculation the advective flux
NEW_PROP_TAG(EnableMolecularDiffusion); //! specifies if molecular diffusive fluxes are considered in the model
NEW_PROP_TAG(MolecularDiffusionType); //! The type for the calculation of the molecular diffusion fluxes
NEW_PROP_TAG(EnableEnergyBalance); //! Specifies if the model solves an energy equation
NEW_PROP_TAG(HeatConductionType); //! The type for the calculation of the heat conduction fluxes
// stencils
......
......@@ -37,10 +37,12 @@
#include <dumux/porousmediumflow/implicit/fluxvariables.hh>
#include <dumux/porousmediumflow/implicit/fluxvariablescache.hh>
#include <dumux/porousmediumflow/nonisothermal/implicit/localresidual.hh>
#include <dumux/discretization/volumevariables.hh>
#include <dumux/discretization/darcyslaw.hh>
#include <dumux/discretization/fickslaw.hh>
#include <dumux/discretization/fourierslaw.hh>
#include "properties.hh"
#include "model.hh"
......@@ -103,20 +105,20 @@ private:
static constexpr bool diffusion = GET_PROP_VALUE(TypeTag, EnableMolecularDiffusion);
static constexpr bool energy = GET_PROP_VALUE(TypeTag, EnableEnergyBalance);
public:
typedef Dumux::PorousMediumFluxVariables<TypeTag, advection, diffusion, energy> type;
typedef PorousMediumFluxVariables<TypeTag, advection, diffusion, energy> type;
};
//! The flux variables cache class, by default the one for porous media
SET_TYPE_PROP(ImplicitBase, FluxVariablesCache, Dumux::PorousMediumFluxVariablesCache<TypeTag>);
SET_TYPE_PROP(ImplicitBase, FluxVariablesCache, PorousMediumFluxVariablesCache<TypeTag>);
//! We use darcys law as the default for the advective flux calculation
SET_TYPE_PROP(ImplicitBase, AdvectionType, Dumux::DarcysLaw<TypeTag>);
//! We use darcy's law as the default for the advective fluxes
SET_TYPE_PROP(ImplicitBase, AdvectionType, DarcysLaw<TypeTag>);
//! We use darcys law as the default for the deffusive flux calculation
SET_TYPE_PROP(ImplicitBase, MolecularDiffusionType, Dumux::FicksLaw<TypeTag>);
//! We use fick's law as the default for the diffusive fluxes
SET_TYPE_PROP(ImplicitBase, MolecularDiffusionType, FicksLaw<TypeTag>);
//! TODO: IMPLEMENT The energy flux çalculation law
//SET_TYPE_PROP(ImplicitBase, HeatConductionType, FouriersLaw<TypeTag>);
//! We use fourier's law as the default for heat conduction fluxes
SET_TYPE_PROP(ImplicitBase, HeatConductionType, FouriersLaw<TypeTag>);
//! The type of a solution for the whole grid at a fixed time
SET_TYPE_PROP(ImplicitBase,
......@@ -164,11 +166,11 @@ SET_BOOL_PROP(ImplicitBase, ConstantBoundaryConditions, false);
SET_PROP(ImplicitBase, JacobianMatrix)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
typedef typename Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
using MatrixBlock = typename Dune::FieldMatrix<Scalar, numEq, numEq>;
public:
typedef typename Dune::BCRSMatrix<MatrixBlock> type;
using type = typename Dune::BCRSMatrix<MatrixBlock>;
};
//! use the stabilized BiCG solver preconditioned by the ILU-0 by default
......@@ -190,6 +192,12 @@ SET_INT_PROP(ImplicitBase, LinearSolverBlockSize, GET_PROP_VALUE(TypeTag, NumEq)
//! set number of maximum timestep divisions to 10
SET_INT_PROP(ImplicitBase, ImplicitMaxTimeStepDivisions, 10);
//! Per default we have assume isothermal problems. Set this to true to solve an energy equation
SET_BOOL_PROP(ImplicitBase, EnableEnergyBalance, false);
SET_TYPE_PROP(ImplicitBase, EnergyLocalResidual, EnergyLocalResidual<TypeTag> );
} // namespace Properties
} // namespace Dumux
......
......@@ -41,23 +41,21 @@ namespace Dumux
template <class TypeTag>
class OnePVolumeVariables : public ImplicitVolumeVariables<TypeTag>
{
typedef ImplicitVolumeVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
using ParentType = ImplicitVolumeVariables<TypeTag>;
using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
public:
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
/*!
* \copydoc ImplicitVolumeVariables::update
......@@ -72,9 +70,6 @@ public:
completeFluidState(priVars, problem, element, scv, fluidState_);
// porosity
porosity_ = problem.spatialParams().porosity(scv);
// energy related quantities not contained in the fluid state
asImp_().updateEnergy_(priVars, problem, element, scv);
};
/*!
......@@ -86,7 +81,7 @@ public:
const SubControlVolume& scv,
FluidState& fluidState)
{
Scalar t = Implementation::temperature_(priVars, problem, element, scv);
Scalar t = ParentType::temperature(priVars, problem, element, scv);
fluidState.setTemperature(t);
fluidState.setSaturation(/*phaseIdx=*/0, 1.);
......@@ -107,7 +102,7 @@ public:
fluidState.setViscosity(/*phaseIdx=*/0, value);
// compute and set the enthalpy
value = Implementation::enthalpy_(fluidState, paramCache, /*phaseIdx=*/0);
value = ParentType::enthalpy(fluidState, paramCache, /*phaseIdx=*/0);
fluidState.setEnthalpy(/*phaseIdx=*/0, value);
}
......@@ -173,35 +168,9 @@ public:
{ return fluidState_; }
protected:
static Scalar temperature_(const PrimaryVariables &priVars,
const Problem& problem,
const Element &element,
const SubControlVolume &scv)
{
return problem.temperatureAtPos(scv.dofPosition());
}
template<class ParameterCache>
static Scalar enthalpy_(const FluidState& fluidState,
const ParameterCache& paramCache,
const int phaseIdx)
{
return 0;
}
/*!
* \brief Called by update() to compute the energy related quantities.
*/
void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const SubControlVolume& scv)
{ }
FluidState fluidState_;
Scalar porosity_;
private:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
......
......@@ -49,7 +49,7 @@ class ImmiscibleLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
// first index for the mass balance
enum { conti0EqIdx = Indices::conti0EqIdx };
......@@ -90,8 +90,14 @@ public:
storage[eqIdx] = volVars.porosity()
* volVars.density(phaseIdx)
* volVars.saturation(phaseIdx);
//! The energy storage in the fluid phase with index phaseIdx
EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
}
//! The energy storage in the solid matrix
EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
return storage;
}
......@@ -127,7 +133,14 @@ public:
auto eqIdx = conti0EqIdx + phaseIdx;
flux[eqIdx] = fluxVars.advectiveFlux(phaseIdx, upwindRule);
//! Add advective phase energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx, w);
}
//! Add diffusive energy fluxes. For isothermal model the contribution is zero.
EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
return flux;
}
......
......@@ -31,8 +31,8 @@ namespace Dumux
namespace Properties
{
// forward declaration
NEW_PROP_TAG(NumPhases);
NEW_PROP_TAG(NumComponents);
}
/*!
......@@ -47,7 +47,8 @@ class PorousMediumFluxVariables {};
// specialization for pure advective flow (e.g. 1p/2p/3p immiscible darcy flow)
template<class TypeTag>
class PorousMediumFluxVariables<TypeTag, true, false, false> : public FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, false, false>>
class PorousMediumFluxVariables<TypeTag, true, false, false>
: public FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, false, false>>
{
using ParentType = FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, false, false>>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
......@@ -106,7 +107,8 @@ public:
// specialization for isothermal advection molecularDiffusion equations
template<class TypeTag>
class PorousMediumFluxVariables<TypeTag, true, true, false> : public FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, true, false>>
class PorousMediumFluxVariables<TypeTag, true, true, false>
: public FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, true, false>>
{
using ParentType = FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, true, false>>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
......@@ -123,11 +125,7 @@ class PorousMediumFluxVariables<TypeTag, true, true, false> : public FluxVariabl
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using MolecularDiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
enum
{
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
};
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
public:
......@@ -206,19 +204,92 @@ private:
// specialization for non-isothermal advective flow (e.g. non-isothermal one-phase darcy equation)
template<class TypeTag>
class PorousMediumFluxVariables<TypeTag, true, false, true>
: public FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, false, true>>
{
using ParentType = FluxVariablesBase<TypeTag, PorousMediumFluxVariables<TypeTag, true, false, true>>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Stencil = std::vector<IndexType>;
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
public:
enum
void initAndComputeFluxes(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& elemVolVars,
const SubControlVolumeFace &scvFace,
const FluxVariablesCache& fluxVarsCache)
{
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
};
advFluxCached_.reset();
ParentType::init(problem, element, fvGeometry, elemVolVars, scvFace, fluxVarsCache);
}
public:
template<typename FunctionType>
Scalar advectiveFlux(const int phaseIdx, const FunctionType& upwindFunction)
{
if (!advFluxCached_[phaseIdx])
{
advPreFlux_[phaseIdx] = AdvectionType::flux(this->problem(),
this->element(),
this->fvGeometry(),
this->elemVolVars(),
this->scvFace(),
phaseIdx,
this->fluxVarsCache());
advFluxCached_.set(phaseIdx, true);
}
const auto& insideScv = this->fvGeometry().scv(this->scvFace().insideScvIdx());
const auto& insideVolVars = this->elemVolVars()[insideScv];
const auto& outsideVolVars = this->elemVolVars()[this->scvFace().outsideScvIdx()];
if (std::signbit(advPreFlux_[phaseIdx]))
return advPreFlux_[phaseIdx]*upwindFunction(outsideVolVars, insideVolVars);
else
return advPreFlux_[phaseIdx]*upwindFunction(insideVolVars, outsideVolVars);
}
Scalar heatConductionFlux()
{
Scalar flux = HeatConductionType::flux(this->problem(),
this->element(),
this->fvGeometry(),
this->elemVolVars(),
this->scvFace(),
this->fluxVarsCache());
return flux;
}
Stencil computeStencil(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvFace)
{
// unifiy advective and diffusive stencil
Stencil stencil = AdvectionType::stencil(problem, element, fvGeometry, scvFace);
Stencil energyStencil = HeatConductionType::stencil(problem, element, fvGeometry, scvFace);
stencil.insert(stencil.end(), energyStencil.begin(), energyStencil.end());
std::sort(stencil.begin(), stencil.end());
stencil.erase(std::unique(stencil.begin(), stencil.end()), stencil.end());
return stencil;
}
private:
//! simple caching if advection flux is used twice with different upwind function
......@@ -229,20 +300,107 @@ private: