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>
151 template <
class LhsEval>
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()) {
247 calculateFluxes_(
flux,
269 RateVector
darcy = 0.0;
271 const auto& problem =
elemCtx.problem();
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,
384 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
387 unsigned pvtRegionIdx =
up.pvtRegionIndex();
395 if constexpr (enableEnergy) {
404 if constexpr (enableEnergy) {
405 EnergyModule::template
414 static_assert(!enableSolvent,
"Relevant computeFlux() method must be implemented for this module before enabling.");
418 static_assert(!enableExtbo,
"Relevant computeFlux() method must be implemented for this module before enabling.");
422 static_assert(!enablePolymer,
"Relevant computeFlux() method must be implemented for this module before enabling.");
426 if constexpr(enableConvectiveMixing) {
427 ConvectiveMixingModule::addConvectiveMixingFlux(
flux,
435 moduleParams.convectiveMixingModuleParam);
439 if constexpr(enableEnergy){
440 const Scalar inAlpha =
nbInfo.inAlpha;
441 const Scalar outAlpha =
nbInfo.outAlpha;
447 EnergyModule::ExtensiveQuantities::template updateEnergy(
heatFlux,
465 static_assert(!enableFoam,
"Relevant computeFlux() method must be implemented for this module before enabling.");
469 static_assert(!enableBrine,
"Relevant computeFlux() method must be implemented for this module before enabling.");
475 if constexpr(enableDiffusion){
476 typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
478 const Scalar diffusivity =
nbInfo.diffusivity;
480 DiffusionModule::addDiffusiveFlux(
flux,
484 effectiveDiffusionCoefficient);
488 if constexpr(enableDispersion){
489 typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg;
491 const Scalar dispersivity =
nbInfo.dispersivity;
493 DispersionModule::addDispersiveFlux(
flux,
501 static_assert(!enableMICP,
"Relevant computeFlux() method must be implemented for this module before enabling.");
506 template <
class BoundaryConditionData>
507 static void computeBoundaryFlux(RateVector&
bdyFlux,
508 const Problem& problem,
509 const BoundaryConditionData&
bdyInfo,
513 if (
bdyInfo.type == BCType::NONE) {
515 }
else if (
bdyInfo.type == BCType::RATE) {
517 }
else if (
bdyInfo.type == BCType::FREE ||
bdyInfo.type == BCType::DIRICHLET) {
519 }
else if (
bdyInfo.type == BCType::THERMAL) {
522 throw std::logic_error(
"Unknown boundary condition type " + std::to_string(
static_cast<int>(
bdyInfo.type)) +
" in computeBoundaryFlux()." );
526 template <
class BoundaryConditionData>
527 static void computeBoundaryFluxRate(RateVector&
bdyFlux,
528 const BoundaryConditionData&
bdyInfo)
533 template <
class BoundaryConditionData>
534 static void computeBoundaryFluxFree(
const Problem& problem,
536 const BoundaryConditionData&
bdyInfo,
541 std::array<short, numPhases>
upIdx;
542 std::array<short, numPhases>
dnIdx;
543 std::array<Evaluation, numPhases> volumeFlux;
544 std::array<Evaluation, numPhases> pressureDifference;
546 ExtensiveQuantities::calculateBoundaryGradients_(problem,
563 if (!FluidSystem::phaseIsActive(
phaseIdx)) {
582 if constexpr (enableEnergy) {
583 EnergyModule::template
589 using ScalarFluidState =
decltype(
bdyInfo.exFluidState);
597 if constexpr (enableEnergy) {
598 EnergyModule::template
607 for (
unsigned i = 0; i < tmp.size(); ++i) {
613 if constexpr(enableEnergy){
616 Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(
globalSpaceIdx,
bdyInfo.boundaryFaceIndex);
619 EnergyModule::ExtensiveQuantities::template updateEnergyBoundary(
heatFlux,
628 static_assert(!enableSolvent,
"Relevant treatment of boundary conditions must be implemented before enabling.");
629 static_assert(!enablePolymer,
"Relevant treatment of boundary conditions must be implemented before enabling.");
630 static_assert(!enableMICP,
"Relevant treatment of boundary conditions must be implemented before enabling.");
636 for (
unsigned i = 0; i < numEq; ++i) {
637 Valgrind::CheckDefined(
bdyFlux[i]);
639 Valgrind::CheckDefined(
bdyFlux);
643 template <
class BoundaryConditionData>
644 static void computeBoundaryThermal(
const Problem& problem,
646 const BoundaryConditionData&
bdyInfo,
655 if constexpr(enableEnergy){
658 Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(
globalSpaceIdx,
bdyInfo.boundaryFaceIndex);
661 EnergyModule::ExtensiveQuantities::template updateEnergyBoundary(
heatFlux,
671 for (
unsigned i = 0; i < numEq; ++i) {
672 Valgrind::CheckDefined(
bdyFlux[i]);
674 Valgrind::CheckDefined(
bdyFlux);
678 static void computeSource(RateVector& source,
679 const Problem& problem,
689 static_assert(!enableMICP,
"Relevant addSource() method must be implemented for this module before enabling.");
697 static void computeSourceDense(RateVector& source,
698 const Problem& problem,
707 static_assert(!enableMICP,
"Relevant addSource() method must be implemented for this module before enabling.");
731 if constexpr(enableEnergy)
735 template <
class UpEval,
class Flu
idState>
736 static void evalPhaseFluxes_(RateVector&
flux,
738 unsigned pvtRegionIdx,
740 const FluidState&
upFs)
752 template <
class UpEval,
class Eval,
class Flu
idState>
755 unsigned pvtRegionIdx,
757 const FluidState&
upFs)
759 unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(
phaseIdx));
761 if (blackoilConserveSurfaceVolume)
768 if (FluidSystem::enableDissolvedGas()) {
769 const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
771 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
772 if (blackoilConserveSurfaceVolume)
777 }
else if (
phaseIdx == waterPhaseIdx) {
779 if (FluidSystem::enableDissolvedGasInWater()) {
780 const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
782 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
783 if (blackoilConserveSurfaceVolume)
791 if (FluidSystem::enableVaporizedOil()) {
792 const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
794 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
795 if (blackoilConserveSurfaceVolume)
801 if (FluidSystem::enableVaporizedWater()) {
802 const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(
upFs, pvtRegionIdx);
805 if (blackoilConserveSurfaceVolume)
824 template <
class Scalar>
827 if (blackoilConserveSurfaceVolume)
837 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
841 unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
843 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
847 unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
849 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
854 static FaceDir::DirEnum faceDirFromDirId(
const int dirId)
858 return FaceDir::DirEnum::Unknown;
860 return FaceDir::FromIntersectionIndex(dirId);
static void computeFlux(RateVector &flux, RateVector &darcy, const unsigned globalIndexIn, const unsigned globalIndexEx, const IntensiveQuantities &intQuantsIn, const IntensiveQuantities &intQuantsEx, const ResidualNBInfo &nbInfo, const ModuleParams &moduleParams)
This function works like the ElementContext-based version with one main difference: The darcy flux is...
Definition blackoillocalresidualtpfa.hh:234