66 using FluidState =
typename IntensiveQuantities::FluidState;
68 enum { conti0EqIdx = Indices::conti0EqIdx };
73 enum { dimWorld = GridView::dimensionworld };
74 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
75 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
76 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
78 enum { gasCompIdx = FluidSystem::gasCompIdx };
79 enum { oilCompIdx = FluidSystem::oilCompIdx };
80 enum { waterCompIdx = FluidSystem::waterCompIdx };
81 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
83 static const bool waterEnabled = Indices::waterEnabled;
84 static const bool gasEnabled = Indices::gasEnabled;
85 static const bool oilEnabled = Indices::oilEnabled;
86 static const bool compositionSwitchEnabled = (compositionSwitchIdx >= 0);
109 using ConvectiveMixingModuleParam =
typename ConvectiveMixingModule::ConvectiveMixingModuleParam;
124 FaceDir::DirEnum faceDir;
134 ConvectiveMixingModuleParam convectiveMixingModuleParam;
140 template <
class LhsEval>
142 const ElementContext& elemCtx,
146 const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx,
timeIdx);
151 template <
class LhsEval>
153 const IntensiveQuantities& intQuants)
157 const auto& fs = intQuants.fluidState();
161 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
164 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
173 if (
phaseIdx == oilPhaseIdx && FluidSystem::enableDissolvedGas()) {
174 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
181 if (
phaseIdx == waterPhaseIdx && FluidSystem::enableDissolvedGasInWater()) {
182 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
189 if (
phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedOil()) {
190 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
197 if (
phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedWater()) {
208 SolventModule::addStorage(
storage, intQuants);
211 ExtboModule::addStorage(
storage, intQuants);
214 PolymerModule::addStorage(
storage, intQuants);
217 EnergyModule::addStorage(
storage, intQuants);
220 FoamModule::addStorage(
storage, intQuants);
223 BrineModule::addStorage(
storage, intQuants);
226 MICPModule::addStorage(
storage, intQuants);
247 calculateFluxes_(
flux,
261 const ElementContext& elemCtx,
269 RateVector
darcy = 0.0;
271 const auto& problem = elemCtx.problem();
272 const auto& stencil = elemCtx.stencil(
timeIdx);
286 Scalar faceArea =
scvf.area();
287 const auto faceDir = faceDirFromDirId(
scvf.dirId());
293 const Scalar
g = problem.gravity()[dimWorld - 1];
313 const ResidualNBInfo res_nbinfo {trans, faceArea, thpres,
distZ *
g, faceDir, Vin, Vex, inAlpha, outAlpha, diffusivity, dispersivity};
315 calculateFluxes_(
flux,
322 problem.moduleParams());
325 static void calculateFluxes_(RateVector&
flux,
331 const ResidualNBInfo&
nbInfo,
332 const ModuleParams& moduleParams)
335 const Scalar Vin =
nbInfo.Vin;
336 const Scalar Vex =
nbInfo.Vex;
338 const Scalar thpres =
nbInfo.thpres;
339 const Scalar trans =
nbInfo.trans;
340 const Scalar faceArea =
nbInfo.faceArea;
344 if (!FluidSystem::phaseIsActive(
phaseIdx))
353 Evaluation pressureDifference;
354 ExtensiveQuantities::calculatePhasePressureDiff_(
upIdx,
376 if constexpr (enableMICP) {
387 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
390 unsigned pvtRegionIdx =
up.pvtRegionIndex();
398 if constexpr (enableEnergy) {
402 if constexpr (enableMICP) {
411 if constexpr (enableEnergy) {
412 EnergyModule::template
416 if constexpr (enableMICP) {
426 static_assert(!enableSolvent,
"Relevant computeFlux() method must be implemented for this module before enabling.");
430 static_assert(!enableExtbo,
"Relevant computeFlux() method must be implemented for this module before enabling.");
434 static_assert(!enablePolymer,
"Relevant computeFlux() method must be implemented for this module before enabling.");
438 if constexpr(enableConvectiveMixing) {
439 ConvectiveMixingModule::addConvectiveMixingFlux(
flux,
447 moduleParams.convectiveMixingModuleParam);
451 if constexpr(enableEnergy){
452 const Scalar inAlpha =
nbInfo.inAlpha;
453 const Scalar outAlpha =
nbInfo.outAlpha;
459 EnergyModule::ExtensiveQuantities::updateEnergy(
heatFlux,
477 static_assert(!enableFoam,
"Relevant computeFlux() method must be implemented for this module before enabling.");
481 static_assert(!enableBrine,
"Relevant computeFlux() method must be implemented for this module before enabling.");
485 if constexpr(enableDiffusion){
486 typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
488 const Scalar diffusivity =
nbInfo.diffusivity;
490 DiffusionModule::addDiffusiveFlux(
flux,
494 effectiveDiffusionCoefficient);
499 if constexpr(enableDispersion){
500 typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg;
502 const Scalar dispersivity =
nbInfo.dispersivity;
504 DispersionModule::addDispersiveFlux(
flux,
513 if constexpr (enableMICP) {
514 MICPModule::applyScaling(
flux);
518 template <
class BoundaryConditionData>
519 static void computeBoundaryFlux(RateVector&
bdyFlux,
520 const Problem& problem,
521 const BoundaryConditionData&
bdyInfo,
525 if (
bdyInfo.type == BCType::NONE) {
527 }
else if (
bdyInfo.type == BCType::RATE) {
529 }
else if (
bdyInfo.type == BCType::FREE ||
bdyInfo.type == BCType::DIRICHLET) {
531 }
else if (
bdyInfo.type == BCType::THERMAL) {
534 throw std::logic_error(
"Unknown boundary condition type " + std::to_string(
static_cast<int>(
bdyInfo.type)) +
" in computeBoundaryFlux()." );
538 template <
class BoundaryConditionData>
539 static void computeBoundaryFluxRate(RateVector&
bdyFlux,
540 const BoundaryConditionData&
bdyInfo)
545 template <
class BoundaryConditionData>
546 static void computeBoundaryFluxFree(
const Problem& problem,
548 const BoundaryConditionData&
bdyInfo,
553 std::array<short, numPhases>
upIdx;
554 std::array<short, numPhases>
dnIdx;
555 std::array<Evaluation, numPhases> volumeFlux;
556 std::array<Evaluation, numPhases> pressureDifference;
558 ExtensiveQuantities::calculateBoundaryGradients_(problem,
575 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
594 if constexpr (enableEnergy) {
595 EnergyModule::template
601 using ScalarFluidState =
decltype(
bdyInfo.exFluidState);
609 if constexpr (enableEnergy) {
610 EnergyModule::template
619 for (
unsigned i = 0; i < tmp.size(); ++i) {
625 if constexpr(enableEnergy){
628 Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(
globalSpaceIdx,
bdyInfo.boundaryFaceIndex);
631 EnergyModule::ExtensiveQuantities::updateEnergyBoundary(
heatFlux,
640 static_assert(!enableSolvent,
"Relevant treatment of boundary conditions must be implemented before enabling.");
641 static_assert(!enablePolymer,
"Relevant treatment of boundary conditions must be implemented before enabling.");
647 for (
unsigned i = 0; i < numEq; ++i) {
648 Valgrind::CheckDefined(
bdyFlux[i]);
650 Valgrind::CheckDefined(
bdyFlux);
654 template <
class BoundaryConditionData>
655 static void computeBoundaryThermal(
const Problem& problem,
657 const BoundaryConditionData&
bdyInfo,
666 if constexpr(enableEnergy){
669 Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(
globalSpaceIdx,
bdyInfo.boundaryFaceIndex);
672 EnergyModule::ExtensiveQuantities::updateEnergyBoundary(
heatFlux,
682 for (
unsigned i = 0; i < numEq; ++i) {
683 Valgrind::CheckDefined(
bdyFlux[i]);
685 Valgrind::CheckDefined(
bdyFlux);
689 static void computeSource(RateVector& source,
690 const Problem& problem,
703 if constexpr(enableEnergy)
707 static void computeSourceDense(RateVector& source,
708 const Problem& problem,
720 if constexpr(enableEnergy)
728 const ElementContext& elemCtx,
734 elemCtx.problem().source(source, elemCtx, dofIdx,
timeIdx);
737 MICPModule::addSource(source, elemCtx, dofIdx,
timeIdx);
740 if constexpr(enableEnergy)
744 template <
class UpEval,
class Flu
idState>
745 static void evalPhaseFluxes_(RateVector&
flux,
747 unsigned pvtRegionIdx,
749 const FluidState&
upFs)
761 template <
class UpEval,
class Eval,
class Flu
idState>
764 unsigned pvtRegionIdx,
766 const FluidState&
upFs)
768 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
770 if (blackoilConserveSurfaceVolume)
777 if (FluidSystem::enableDissolvedGas()) {
778 const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
780 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
781 if (blackoilConserveSurfaceVolume)
786 }
else if (
phaseIdx == waterPhaseIdx) {
788 if (FluidSystem::enableDissolvedGasInWater()) {
789 const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
791 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
792 if (blackoilConserveSurfaceVolume)
800 if (FluidSystem::enableVaporizedOil()) {
801 const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
803 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
804 if (blackoilConserveSurfaceVolume)
810 if (FluidSystem::enableVaporizedWater()) {
811 const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
814 if (blackoilConserveSurfaceVolume)
833 template <
class Scalar>
836 if (blackoilConserveSurfaceVolume)
846 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
850 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
852 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
856 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
858 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
863 static FaceDir::DirEnum faceDirFromDirId(
const int dirId)
867 return FaceDir::DirEnum::Unknown;
869 return FaceDir::FromIntersectionIndex(dirId);