28#ifndef EWOMS_BLACK_OIL_ENERGY_MODULE_HH
29#define EWOMS_BLACK_OIL_ENERGY_MODULE_HH
31#include <dune/common/fvector.hh>
33#include <opm/material/common/Tabulated1DFunction.hpp>
34#include <opm/material/common/Valgrind.hpp>
56template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()>
71 static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
72 static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;
74 static constexpr unsigned enableEnergy = enableEnergyV;
75 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
76 static constexpr unsigned numPhases = FluidSystem::numPhases;
86 if constexpr (enableEnergy) {
97 if constexpr (enableEnergy) {
104 if constexpr (enableEnergy) {
105 return pvIdx == temperatureIdx;
116 return "temperature";
124 return static_cast<Scalar
>(1.0);
129 if constexpr (enableEnergy) {
130 return eqIdx == contiEnergyEqIdx;
137 static std::string
eqName([[maybe_unused]]
unsigned eqIdx)
141 return "conti^energy";
144 static Scalar
eqWeight([[maybe_unused]]
unsigned eqIdx)
152 template <
class LhsEval>
153 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
154 const IntensiveQuantities& intQuants)
156 if constexpr (enableEnergy) {
157 const auto& poro = decay<LhsEval>(intQuants.porosity());
160 const auto& fs = intQuants.fluidState();
161 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
162 if (!FluidSystem::phaseIsActive(phaseIdx)) {
166 const auto& u = decay<LhsEval>(fs.internalEnergy(phaseIdx));
167 const auto& S = decay<LhsEval>(fs.saturation(phaseIdx));
168 const auto& rho = decay<LhsEval>(fs.density(phaseIdx));
170 storage[contiEnergyEqIdx] += poro*S*u*rho;
174 const Scalar rockFraction = intQuants.rockFraction();
175 const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
176 storage[contiEnergyEqIdx] += rockFraction * uRock;
177 storage[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
182 [[maybe_unused]]
const ElementContext& elemCtx,
183 [[maybe_unused]]
unsigned scvfIdx,
184 [[maybe_unused]]
unsigned timeIdx)
186 if constexpr (enableEnergy) {
187 flux[contiEnergyEqIdx] = 0.0;
189 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
190 const unsigned focusIdx = elemCtx.focusDofIndex();
191 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
192 if (!FluidSystem::phaseIsActive(phaseIdx)) {
196 const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
197 if (upIdx == focusIdx) {
198 addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
201 addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
206 flux[contiEnergyEqIdx] += extQuants.energyFlux();
207 flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
212 const Evaluation& heatFlux)
214 if constexpr (enableEnergy) {
216 flux[contiEnergyEqIdx] += heatFlux;
217 flux[contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
221 template <
class UpEval,
class Eval,
class Flu
idState>
224 const Eval& volumeFlux,
225 const FluidState& upFs)
227 flux[contiEnergyEqIdx] +=
228 decay<UpEval>(upFs.enthalpy(phaseIdx)) *
229 decay<UpEval>(upFs.density(phaseIdx)) *
233 template <
class UpstreamEval>
236 const ElementContext& elemCtx,
240 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
241 const unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
242 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
243 const auto& fs = up.fluidState();
244 const auto& volFlux = extQuants.volumeFlux(phaseIdx);
245 addPhaseEnthalpyFluxes_<UpstreamEval>(flux,
252 const Evaluation& hRate)
254 if constexpr (enableEnergy) {
255 flux[contiEnergyEqIdx] += hRate;
262 template <
class Flu
idState>
264 const FluidState& fluidState)
266 if constexpr (enableEnergy) {
267 priVars[temperatureIdx] = fluidState.temperature(0);
275 const PrimaryVariables& oldPv,
276 const EqVector& delta)
278 if constexpr (enableEnergy) {
280 newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
293 return static_cast<Scalar
>(0.0);
302 return std::abs(scalarValue(resid[contiEnergyEqIdx]));
305 template <
class DofEntity>
306 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
308 if constexpr (enableEnergy) {
309 const unsigned dofIdx = model.dofMapper().index(dof);
310 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
311 outstream << priVars[temperatureIdx];
315 template <
class DofEntity>
318 if constexpr (enableEnergy) {
319 const unsigned dofIdx = model.dofMapper().index(dof);
320 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
321 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
323 instream >> priVars0[temperatureIdx];
326 priVars1 = priVars0[temperatureIdx];
331template <
class TypeTag,
bool enableEnergyV>
341template <
class TypeTag>
356 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
357 static constexpr int temperatureIdx = Indices::temperatureIdx;
368 auto& fs = asImp_().fluidState_;
369 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
372 fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, elemCtx.linearizationType()));
380 const PrimaryVariables& priVars,
381 [[maybe_unused]]
unsigned globalDofIdx,
382 const unsigned timeIdx,
385 auto& fs = asImp_().fluidState_;
386 fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, lintype));
397 updateEnergyQuantities_(elemCtx.problem(), elemCtx.globalSpaceIndex(dofIdx, timeIdx), timeIdx);
401 const unsigned globalSpaceIdx,
402 const unsigned timeIdx)
404 auto& fs = asImp_().fluidState_;
408 for (
int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
409 if (!FluidSystem::phaseIsActive(phaseIdx)) {
413 const auto& h = FluidSystem::enthalpy(fs, phaseIdx, problem.pvtRegionIndex(globalSpaceIdx));
414 fs.setEnthalpy(phaseIdx, h);
417 const auto& solidEnergyLawParams = problem.solidEnergyLawParams(globalSpaceIdx, timeIdx);
418 rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
420 const auto& thermalConductionLawParams = problem.thermalConductionLawParams(globalSpaceIdx, timeIdx);
421 totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
428 rockFraction_ = problem.rockFraction(globalSpaceIdx, timeIdx);
432 {
return rockInternalEnergy_; }
435 {
return totalThermalConductivity_; }
438 {
return rockFraction_; }
442 {
return *
static_cast<Implementation*
>(
this); }
449template <
class TypeTag>
460 static constexpr bool enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>();
464 [[maybe_unused]]
unsigned dofIdx,
465 [[maybe_unused]]
unsigned timeIdx)
467 if constexpr (enableTemperature) {
470 auto& fs = asImp_().fluidState_;
471 const Scalar T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
472 fs.setTemperature(T);
476 template<
class Problem>
478 [[maybe_unused]]
const PrimaryVariables& priVars,
479 [[maybe_unused]]
unsigned globalDofIdx,
480 [[maybe_unused]]
unsigned timeIdx,
484 if constexpr (enableTemperature) {
485 auto& fs = asImp_().fluidState_;
488 const Scalar T = problem.temperature(globalDofIdx, timeIdx);
489 fs.setTemperature(T);
496 const typename FluidSystem::template ParameterCache<Evaluation>&)
501 throw std::logic_error(
"Requested the rock internal energy, which is "
502 "unavailable because energy is not conserved");
507 throw std::logic_error(
"Requested the total thermal conductivity, which is "
508 "unavailable because energy is not conserved");
513 {
return *
static_cast<Implementation*
>(
this); }
516template <
class TypeTag,
bool enableEnergyV>
526template <
class TypeTag>
537 template<
class Flu
idState>
539 const unsigned& focusDofIndex,
540 const unsigned& inIdx,
541 const unsigned& exIdx,
542 const IntensiveQuantities& inIq,
543 const IntensiveQuantities& exIq,
544 const FluidState& inFs,
545 const FluidState& exFs,
546 const Scalar& inAlpha,
547 const Scalar& outAlpha,
548 const Scalar& faceArea)
551 if (focusDofIndex == inIdx) {
552 deltaT = decay<Scalar>(exFs.temperature(0)) -
555 else if (focusDofIndex == exIdx) {
556 deltaT = exFs.temperature(0) -
557 decay<Scalar>(inFs.temperature(0));
560 deltaT = decay<Scalar>(exFs.temperature(0)) -
561 decay<Scalar>(inFs.temperature(0));
565 if (focusDofIndex == inIdx) {
566 inLambda = inIq.totalThermalConductivity();
569 inLambda = decay<Scalar>(inIq.totalThermalConductivity());
573 if (focusDofIndex == exIdx) {
574 exLambda = exIq.totalThermalConductivity();
577 exLambda = decay<Scalar>(exIq.totalThermalConductivity());
581 const Evaluation& inH = inLambda*inAlpha;
582 const Evaluation& exH = exLambda*outAlpha;
583 if (inH > 0 && exH > 0) {
588 H = 1.0 / (1.0 / inH + 1.0 / exH);
594 energyFlux = deltaT * (-H / faceArea);
601 const auto& stencil = elemCtx.stencil(timeIdx);
602 const auto& scvf = stencil.interiorFace(scvfIdx);
604 const Scalar faceArea = scvf.area();
605 const unsigned inIdx = scvf.interiorIndex();
606 const unsigned exIdx = scvf.exteriorIndex();
607 const auto& inIq = elemCtx.intensiveQuantities(inIdx, timeIdx);
608 const auto& exIq = elemCtx.intensiveQuantities(exIdx, timeIdx);
609 const auto& inFs = inIq.fluidState();
610 const auto& exFs = exIq.fluidState();
611 const Scalar inAlpha = elemCtx.problem().thermalHalfTransmissibilityIn(elemCtx, scvfIdx, timeIdx);
612 const Scalar outAlpha = elemCtx.problem().thermalHalfTransmissibilityOut(elemCtx, scvfIdx, timeIdx);
613 updateEnergy(energyFlux_,
614 elemCtx.focusDofIndex(),
626 template <
class Context,
class BoundaryFlu
idState>
630 const BoundaryFluidState& boundaryFs)
632 const auto& stencil = ctx.stencil(timeIdx);
633 const auto& scvf = stencil.boundaryFace(scvfIdx);
635 const unsigned inIdx = scvf.interiorIndex();
636 const auto& inIq = ctx.intensiveQuantities(inIdx, timeIdx);
637 const auto& focusDofIdx = ctx.focusDofIndex();
638 const Scalar alpha = ctx.problem().thermalHalfTransmissibilityBoundary(ctx, scvfIdx);
639 updateEnergyBoundary(energyFlux_, inIq, focusDofIdx, inIdx, alpha, boundaryFs);
642 template <
class BoundaryFlu
idState>
644 const IntensiveQuantities& inIq,
645 unsigned focusDofIndex,
648 const BoundaryFluidState& boundaryFs)
650 const auto& inFs = inIq.fluidState();
652 if (focusDofIndex == inIdx) {
653 deltaT = boundaryFs.temperature(0) -
657 deltaT = decay<Scalar>(boundaryFs.temperature(0)) -
658 decay<Scalar>(inFs.temperature(0));
662 if (focusDofIndex == inIdx) {
663 lambda = inIq.totalThermalConductivity();
666 lambda = decay<Scalar>(inIq.totalThermalConductivity());
674 energyFlux = deltaT * lambda * -alpha;
682 {
return energyFlux_; }
685 Implementation& asImp_()
686 {
return *
static_cast<Implementation*
>(
this); }
688 Evaluation energyFlux_;
691template <
class TypeTag>
700 template<
class Flu
idState>
705 const IntensiveQuantities& ,
706 const IntensiveQuantities& ,
719 template <
class Context,
class BoundaryFlu
idState>
723 const BoundaryFluidState&)
726 template <
class BoundaryFlu
idState>
728 const IntensiveQuantities& ,
733 const BoundaryFluidState& )
737 {
throw std::logic_error(
"Requested the energy flux, but energy is not conserved"); }
Declares the properties required by the black oil model.
void updateEnergyBoundary(const Context &, unsigned, unsigned, const BoundaryFluidState &)
Definition: blackoilenergymodules.hh:720
static void updateEnergy(Evaluation &, const unsigned &, const unsigned &, const unsigned &, const IntensiveQuantities &, const IntensiveQuantities &, const FluidState &, const FluidState &, const Scalar &, const Scalar &, const Scalar &)
Definition: blackoilenergymodules.hh:701
static void updateEnergyBoundary(Evaluation &, const IntensiveQuantities &, unsigned, unsigned, unsigned, Scalar, const BoundaryFluidState &)
Definition: blackoilenergymodules.hh:727
void updateEnergy(const ElementContext &, unsigned, unsigned)
Definition: blackoilenergymodules.hh:714
const Evaluation & energyFlux() const
Definition: blackoilenergymodules.hh:736
void updateEnergyBoundary(const Context &ctx, unsigned scvfIdx, unsigned timeIdx, const BoundaryFluidState &boundaryFs)
Definition: blackoilenergymodules.hh:627
static void updateEnergyBoundary(Evaluation &energyFlux, const IntensiveQuantities &inIq, unsigned focusDofIndex, unsigned inIdx, Scalar alpha, const BoundaryFluidState &boundaryFs)
Definition: blackoilenergymodules.hh:643
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)
Definition: blackoilenergymodules.hh:538
const Evaluation & energyFlux() const
Definition: blackoilenergymodules.hh:681
void updateEnergy(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilenergymodules.hh:597
Provides the energy specific extensive quantities to the generic black-oil module's extensive quantit...
Definition: blackoilenergymodules.hh:517
void updateEnergyQuantities_(const ElementContext &, unsigned, unsigned, const typename FluidSystem::template ParameterCache< Evaluation > &)
Definition: blackoilenergymodules.hh:493
const Evaluation & rockInternalEnergy() const
Definition: blackoilenergymodules.hh:499
const Evaluation & totalThermalConductivity() const
Definition: blackoilenergymodules.hh:505
void updateTemperature_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: blackoilenergymodules.hh:463
void updateTemperature_(const Problem &problem, const PrimaryVariables &priVars, unsigned globalDofIdx, unsigned timeIdx, const LinearizationType &lintype)
Definition: blackoilenergymodules.hh:477
Implementation & asImp_()
Definition: blackoilenergymodules.hh:512
Evaluation totalThermalConductivity_
Definition: blackoilenergymodules.hh:445
const Evaluation & totalThermalConductivity() const
Definition: blackoilenergymodules.hh:434
void updateEnergyQuantities_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Compute the intensive quantities needed to handle energy conservation.
Definition: blackoilenergymodules.hh:393
Scalar rockFraction() const
Definition: blackoilenergymodules.hh:437
Implementation & asImp_()
Definition: blackoilenergymodules.hh:441
void updateTemperature_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the temperature of the intensive quantity's fluid state.
Definition: blackoilenergymodules.hh:364
Scalar rockFraction_
Definition: blackoilenergymodules.hh:446
void updateTemperature_(const Problem &problem, const PrimaryVariables &priVars, unsigned globalDofIdx, const unsigned timeIdx, const LinearizationType &lintype)
Update the temperature of the intensive quantity's fluid state.
Definition: blackoilenergymodules.hh:379
const Evaluation & rockInternalEnergy() const
Definition: blackoilenergymodules.hh:431
void updateEnergyQuantities_(const Problem &problem, const unsigned globalSpaceIdx, const unsigned timeIdx)
Definition: blackoilenergymodules.hh:400
Evaluation rockInternalEnergy_
Definition: blackoilenergymodules.hh:444
Provides the volumetric quantities required for the equations needed by the energys extension of the ...
Definition: blackoilenergymodules.hh:332
Contains the high level supplements required to extend the black oil model by energy.
Definition: blackoilenergymodules.hh:58
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilenergymodules.hh:181
static void registerParameters()
Register all run-time parameters for the black-oil energy module.
Definition: blackoilenergymodules.hh:84
GetPropType< TypeTag, Properties::ExtensiveQuantities > ExtensiveQuantities
Definition: blackoilenergymodules.hh:79
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilenergymodules.hh:299
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilenergymodules.hh:316
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilenergymodules.hh:153
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilenergymodules.hh:102
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilenergymodules.hh:112
static std::string eqName(unsigned eqIdx)
Definition: blackoilenergymodules.hh:137
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilenergymodules.hh:306
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilenergymodules.hh:119
static void updatePrimaryVars(PrimaryVariables &newPv, const PrimaryVariables &oldPv, const EqVector &delta)
Do a Newton-Raphson update the primary variables of the energys.
Definition: blackoilenergymodules.hh:274
static void addToEnthalpyRate(RateVector &flux, const Evaluation &hRate)
Definition: blackoilenergymodules.hh:251
static void addHeatFlux(RateVector &flux, const Evaluation &heatFlux)
Definition: blackoilenergymodules.hh:211
static void registerOutputModules(Model &model, Simulator &simulator)
Register all energy specific VTK and ECL output modules.
Definition: blackoilenergymodules.hh:94
static void addPhaseEnthalpyFlux_(RateVector &flux, unsigned phaseIdx, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilenergymodules.hh:234
static void addPhaseEnthalpyFluxes_(RateVector &flux, unsigned phaseIdx, const Eval &volumeFlux, const FluidState &upFs)
Definition: blackoilenergymodules.hh:222
static bool eqApplies(unsigned eqIdx)
Definition: blackoilenergymodules.hh:127
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilenergymodules.hh:287
static void assignPrimaryVars(PrimaryVariables &priVars, const FluidState &fluidState)
Assign the energy specific primary variables to a PrimaryVariables object.
Definition: blackoilenergymodules.hh:263
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilenergymodules.hh:144
VTK output module for the black oil model's energy related quantities.
Definition: vtkblackoilenergymodule.hpp:54
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkblackoilenergymodule.hpp:87
Definition: blackoilbioeffectsmodules.hh:43
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:233
This method contains all callback classes for quantities that are required by some extensive quantiti...
Definition: linearizationtype.hh:34