28#ifndef EWOMS_BLACK_OIL_ENERGY_MODULE_HH
29#define EWOMS_BLACK_OIL_ENERGY_MODULE_HH
36#include <opm/material/common/Tabulated1DFunction.hpp>
38#include <opm/material/common/Valgrind.hpp>
40#include <dune/common/fvector.hh>
50template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()>
65 static constexpr unsigned temperatureIdx = Indices::temperatureIdx;
66 static constexpr unsigned contiEnergyEqIdx = Indices::contiEnergyEqIdx;
68 static constexpr unsigned enableEnergy = enableEnergyV;
70 static constexpr unsigned numPhases = FluidSystem::numPhases;
79 if constexpr (enableEnergy)
89 if constexpr (enableEnergy)
93 static bool primaryVarApplies(
unsigned pvIdx)
95 if constexpr (enableEnergy)
96 return pvIdx == temperatureIdx;
102 static std::string primaryVarName([[maybe_unused]]
unsigned pvIdx)
104 assert(primaryVarApplies(pvIdx));
106 return "temperature";
109 static Scalar primaryVarWeight([[maybe_unused]]
unsigned pvIdx)
111 assert(primaryVarApplies(pvIdx));
114 return static_cast<Scalar
>(1.0);
117 static bool eqApplies(
unsigned eqIdx)
119 if constexpr (enableEnergy)
120 return eqIdx == contiEnergyEqIdx;
125 static std::string eqName([[maybe_unused]]
unsigned eqIdx)
127 assert(eqApplies(eqIdx));
129 return "conti^energy";
132 static Scalar eqWeight([[maybe_unused]]
unsigned eqIdx)
134 assert(eqApplies(eqIdx));
140 template <
class LhsEval>
141 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
142 const IntensiveQuantities& intQuants)
144 if constexpr (enableEnergy) {
145 const auto& poro = decay<LhsEval>(intQuants.porosity());
148 const auto& fs = intQuants.fluidState();
149 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
150 if (!FluidSystem::phaseIsActive(phaseIdx))
153 const auto& u = decay<LhsEval>(fs.internalEnergy(phaseIdx));
154 const auto& S = decay<LhsEval>(fs.saturation(phaseIdx));
155 const auto& rho = decay<LhsEval>(fs.density(phaseIdx));
157 storage[contiEnergyEqIdx] += poro*S*u*rho;
161 Scalar rockFraction = intQuants.rockFraction();
162 const auto& uRock = decay<LhsEval>(intQuants.rockInternalEnergy());
163 storage[contiEnergyEqIdx] += rockFraction*uRock;
168 static void computeFlux([[maybe_unused]] RateVector& flux,
169 [[maybe_unused]]
const ElementContext& elemCtx,
170 [[maybe_unused]]
unsigned scvfIdx,
171 [[maybe_unused]]
unsigned timeIdx)
173 if constexpr (enableEnergy) {
174 flux[contiEnergyEqIdx] = 0.0;
176 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
177 unsigned focusIdx = elemCtx.focusDofIndex();
178 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
179 if (!FluidSystem::phaseIsActive(phaseIdx))
182 unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
183 if (upIdx == focusIdx)
184 addPhaseEnthalpyFlux_<Evaluation>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
186 addPhaseEnthalpyFlux_<Scalar>(flux, phaseIdx, elemCtx, scvfIdx, timeIdx);
190 flux[contiEnergyEqIdx] += extQuants.energyFlux();
195 static void addHeatFlux(RateVector& flux,
196 const Evaluation& heatFlux)
198 if constexpr (enableEnergy) {
200 flux[contiEnergyEqIdx] += heatFlux;
207 template <
class UpEval,
class Eval,
class Flu
idState>
208 static void addPhaseEnthalpyFluxes_(RateVector& flux,
210 const Eval& volumeFlux,
211 const FluidState& upFs)
213 flux[contiEnergyEqIdx] +=
214 decay<UpEval>(upFs.enthalpy(phaseIdx))
215 * decay<UpEval>(upFs.density(phaseIdx))
219 template <
class UpstreamEval>
220 static void addPhaseEnthalpyFlux_(RateVector& flux,
222 const ElementContext& elemCtx,
226 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
227 unsigned upIdx = extQuants.upstreamIndex(phaseIdx);
228 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
229 const auto& fs = up.fluidState();
230 const auto& volFlux = extQuants.volumeFlux(phaseIdx);
231 addPhaseEnthalpyFluxes_<UpstreamEval>(flux,
237 static void addToEnthalpyRate(RateVector& flux,
238 const Evaluation& hRate)
240 if constexpr (enableEnergy)
241 flux[contiEnergyEqIdx] += hRate;
250 if constexpr (enableEnergy)
251 priVars[temperatureIdx] = temperatureIdx;
257 template <
class Flu
idState>
259 const FluidState& fluidState)
261 if constexpr (enableEnergy)
262 priVars[temperatureIdx] = fluidState.temperature(0);
269 const PrimaryVariables& oldPv,
270 const EqVector& delta)
272 if constexpr (enableEnergy)
274 newPv[temperatureIdx] = oldPv[temperatureIdx] - delta[temperatureIdx];
286 return static_cast<Scalar
>(0.0);
295 return std::abs(scalarValue(resid[contiEnergyEqIdx]));
298 template <
class DofEntity>
299 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
301 if constexpr (enableEnergy) {
302 unsigned dofIdx = model.dofMapper().index(dof);
303 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
304 outstream << priVars[temperatureIdx];
308 template <
class DofEntity>
309 static void deserializeEntity(Model& model, std::istream& instream,
const DofEntity& dof)
311 if constexpr (enableEnergy) {
312 unsigned dofIdx = model.dofMapper().index(dof);
313 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
314 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
316 instream >> priVars0[temperatureIdx];
319 priVars1 = priVars0[temperatureIdx];
331template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()>
332class BlackOilEnergyIntensiveQuantities
334 using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
336 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
337 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
338 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
339 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
340 using SolidEnergyLaw = GetPropType<TypeTag, Properties::SolidEnergyLaw>;
341 using ThermalConductionLaw = GetPropType<TypeTag, Properties::ThermalConductionLaw>;
342 using Indices = GetPropType<TypeTag, Properties::Indices>;
343 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
344 using Problem = GetPropType<TypeTag, Properties::Problem>;
346 using EnergyModule = BlackOilEnergyModule<TypeTag>;
349 static constexpr int temperatureIdx = Indices::temperatureIdx;
350 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
362 auto& fs = asImp_().fluidState_;
363 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
366 fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, elemCtx.linearizationType()));
374 const PrimaryVariables& priVars,
375 [[maybe_unused]]
unsigned globalDofIdx,
376 const unsigned timeIdx,
379 auto& fs = asImp_().fluidState_;
380 fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, lintype));
390 const typename FluidSystem::template ParameterCache<Evaluation>& paramCache)
392 auto& fs = asImp_().fluidState_;
396 for (
int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
397 if (!FluidSystem::phaseIsActive(phaseIdx)) {
401 const auto& h = FluidSystem::enthalpy(fs, paramCache, phaseIdx);
402 fs.setEnthalpy(phaseIdx, h);
405 const auto& solidEnergyLawParams = elemCtx.problem().solidEnergyLawParams(elemCtx, dofIdx, timeIdx);
406 rockInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyLawParams, fs);
408 const auto& thermalConductionLawParams = elemCtx.problem().thermalConductionLawParams(elemCtx, dofIdx, timeIdx);
409 totalThermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalConductionLawParams, fs);
416 const unsigned cell_idx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
417 rockFraction_ = elemCtx.problem().rockFraction(cell_idx, timeIdx);
420 const Evaluation& rockInternalEnergy()
const
421 {
return rockInternalEnergy_; }
423 const Evaluation& totalThermalConductivity()
const
424 {
return totalThermalConductivity_; }
426 const Scalar& rockFraction()
const
427 {
return rockFraction_; }
430 Implementation& asImp_()
431 {
return *
static_cast<Implementation*
>(
this); }
433 Evaluation rockInternalEnergy_;
434 Evaluation totalThermalConductivity_;
435 Scalar rockFraction_;
438template <
class TypeTag>
453 [[maybe_unused]]
unsigned dofIdx,
454 [[maybe_unused]]
unsigned timeIdx)
456 if constexpr (enableTemperature) {
459 auto& fs = asImp_().fluidState_;
460 const Scalar T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
461 fs.setTemperature(T);
465 template<
class Problem>
467 [[maybe_unused]]
const PrimaryVariables& priVars,
468 [[maybe_unused]]
unsigned globalDofIdx,
469 [[maybe_unused]]
unsigned timeIdx,
473 if constexpr (enableTemperature) {
474 auto& fs = asImp_().fluidState_;
477 const Scalar T = problem.temperature(globalDofIdx, timeIdx);
478 fs.setTemperature(T);
485 const typename FluidSystem::template ParameterCache<Evaluation>&)
488 const Evaluation& rockInternalEnergy()
const
489 {
throw std::logic_error(
"Requested the rock internal energy, which is "
490 "unavailable because energy is not conserved"); }
492 const Evaluation& totalThermalConductivity()
const
493 {
throw std::logic_error(
"Requested the total thermal conductivity, which is "
494 "unavailable because energy is not conserved"); }
497 Implementation& asImp_()
498 {
return *
static_cast<Implementation*
>(
this); }
509template <class TypeTag, bool enableEnergyV = getPropValue<TypeTag, Properties::EnableEnergy>()>
522 using Toolbox = MathToolbox<Evaluation>;
526 static const int dimWorld = GridView::dimensionworld;
527 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
528 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
530 template<
class Flu
idState>
531 static void updateEnergy(Evaluation& energyFlux,
532 const unsigned& focusDofIndex,
533 const unsigned& inIdx,
534 const unsigned& exIdx,
535 const IntensiveQuantities& inIq,
536 const IntensiveQuantities& exIq,
537 const FluidState& inFs,
538 const FluidState& exFs,
539 const Scalar& inAlpha,
540 const Scalar& outAlpha,
541 const Scalar& faceArea)
544 if (focusDofIndex == inIdx)
546 decay<Scalar>(exFs.temperature(0))
547 - inFs.temperature(0);
548 else if (focusDofIndex == exIdx)
551 - decay<Scalar>(inFs.temperature(0));
554 decay<Scalar>(exFs.temperature(0))
555 - decay<Scalar>(inFs.temperature(0));
558 if (focusDofIndex == inIdx)
559 inLambda = inIq.totalThermalConductivity();
561 inLambda = decay<Scalar>(inIq.totalThermalConductivity());
564 if (focusDofIndex == exIdx)
565 exLambda = exIq.totalThermalConductivity();
567 exLambda = decay<Scalar>(exIq.totalThermalConductivity());
570 const Evaluation& inH = inLambda*inAlpha;
571 const Evaluation& exH = exLambda*outAlpha;
572 if (inH > 0 && exH > 0) {
577 H = 1.0/(1.0/inH + 1.0/exH);
582 energyFlux = deltaT * (-H/faceArea);
585 void updateEnergy(
const ElementContext& elemCtx,
589 const auto& stencil = elemCtx.stencil(timeIdx);
590 const auto& scvf = stencil.interiorFace(scvfIdx);
592 const Scalar faceArea = scvf.area();
593 const unsigned inIdx = scvf.interiorIndex();
594 const unsigned exIdx = scvf.exteriorIndex();
595 const auto& inIq = elemCtx.intensiveQuantities(inIdx, timeIdx);
596 const auto& exIq = elemCtx.intensiveQuantities(exIdx, timeIdx);
597 const auto& inFs = inIq.fluidState();
598 const auto& exFs = exIq.fluidState();
599 const Scalar inAlpha = elemCtx.problem().thermalHalfTransmissibilityIn(elemCtx, scvfIdx, timeIdx);
600 const Scalar outAlpha = elemCtx.problem().thermalHalfTransmissibilityOut(elemCtx, scvfIdx, timeIdx);
601 updateEnergy(energyFlux_,
602 elemCtx.focusDofIndex(),
614 template <
class Context,
class BoundaryFlu
idState>
615 void updateEnergyBoundary(
const Context& ctx,
618 const BoundaryFluidState& boundaryFs)
620 const auto& stencil = ctx.stencil(timeIdx);
621 const auto& scvf = stencil.boundaryFace(scvfIdx);
623 unsigned inIdx = scvf.interiorIndex();
624 const auto& inIq = ctx.intensiveQuantities(inIdx, timeIdx);
625 const auto& focusDofIdx = ctx.focusDofIndex();
626 const Scalar alpha = ctx.problem().thermalHalfTransmissibilityBoundary(ctx, scvfIdx);
627 updateEnergyBoundary(energyFlux_, inIq, focusDofIdx, inIdx, alpha, boundaryFs);
630 template <
class BoundaryFlu
idState>
631 static void updateEnergyBoundary(Evaluation& energyFlux,
632 const IntensiveQuantities& inIq,
633 unsigned focusDofIndex,
636 const BoundaryFluidState& boundaryFs)
638 const auto& inFs = inIq.fluidState();
640 if (focusDofIndex == inIdx)
642 boundaryFs.temperature(0)
643 - inFs.temperature(0);
646 decay<Scalar>(boundaryFs.temperature(0))
647 - decay<Scalar>(inFs.temperature(0));
650 if (focusDofIndex == inIdx)
651 lambda = inIq.totalThermalConductivity();
653 lambda = decay<Scalar>(inIq.totalThermalConductivity());
661 energyFlux = deltaT*lambda*(-alpha);
667 const Evaluation& energyFlux()
const
668 {
return energyFlux_; }
671 Implementation& asImp_()
672 {
return *
static_cast<Implementation*
>(
this); }
674 Evaluation energyFlux_;
677template <
class TypeTag>
685 template<
class Flu
idState>
686 static void updateEnergy(Evaluation& ,
690 const IntensiveQuantities& ,
691 const IntensiveQuantities& ,
699 void updateEnergy(
const ElementContext&,
704 template <
class Context,
class BoundaryFlu
idState>
705 void updateEnergyBoundary(
const Context&,
708 const BoundaryFluidState&)
711 template <
class BoundaryFlu
idState>
712 static void updateEnergyBoundary(Evaluation& ,
713 const IntensiveQuantities& ,
718 const BoundaryFluidState& )
720 const Evaluation& energyFlux()
const
721 {
throw std::logic_error(
"Requested the energy flux, but energy is not conserved"); }
Declares the properties required by the black oil model.
Provides the energy specific extensive quantities to the generic black-oil module's extensive quantit...
Provides the volumetric quantities required for the equations needed by the energys extension of the ...
void updateTemperature_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the temperature of the intensive quantity's fluid state.
Definition blackoilenergymodules.hh:358
void updateEnergyQuantities_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx, const typename FluidSystem::template ParameterCache< Evaluation > ¶mCache)
Compute the intensive quantities needed to handle energy conservation.
Definition blackoilenergymodules.hh:387
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:373
Contains the high level supplements required to extend the black oil model by energy.
Definition blackoilenergymodules.hh:52
static void registerParameters()
Register all run-time parameters for the black-oil energy module.
Definition blackoilenergymodules.hh:77
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition blackoilenergymodules.hh:292
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:268
static void registerOutputModules(Model &model, Simulator &simulator)
Register all energy specific VTK and ECL output modules.
Definition blackoilenergymodules.hh:86
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar)
Assign the energy specific primary variables to a PrimaryVariables object.
Definition blackoilenergymodules.hh:247
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition blackoilenergymodules.hh:280
static void assignPrimaryVars(PrimaryVariables &priVars, const FluidState &fluidState)
Assign the energy specific primary variables to a PrimaryVariables object.
Definition blackoilenergymodules.hh:258
Provides the auxiliary methods required for consideration of the energy equation.
Definition energymodule.hh:50
VTK output module for the black oil model's energy related quantities.
Definition vtkblackoilenergymodule.hpp:54
The common code for the linearizers of non-linear systems of equations.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
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:235
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
This method contains all callback classes for quantities that are required by some extensive quantiti...
Definition linearizationtype.hh:35
VTK output module for the black oil model's energy related quantities.