321 const MaterialLawParams& matParams,
322 bool isInEquilibrium =
false)
324 using ConstEvaluation =
typename std::remove_reference<typename FluidState::Scalar>::type;
325 using FsEvaluation =
typename std::remove_const<ConstEvaluation>::type;
326 using FsToolbox = MathToolbox<FsEvaluation>;
330 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
331 Valgrind::CheckDefined(fluidState.temperature(0));
332 Valgrind::CheckDefined(fluidState.temperature(phaseIdx));
334 assert(fluidState.temperature(0) == fluidState.temperature(phaseIdx));
340 if (isInEquilibrium) {
347 typename FluidSystem::template ParameterCache<Scalar> paramCache;
348 paramCache.setRegionIndex(pvtRegionIdx_);
349 paramCache.setMaxOilSat(FsToolbox::value(fluidState.saturation(oilPhaseIdx)));
352 using NcpFlash = NcpFlash<Scalar, FluidSystem>;
353 using FlashFluidState = CompositionalFluidState<Scalar, FluidSystem>;
354 FlashFluidState fsFlash;
355 fsFlash.setTemperature(FsToolbox::value(fluidState.temperature(0)));
356 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
357 fsFlash.setPressure(phaseIdx, FsToolbox::value(fluidState.pressure(phaseIdx)));
358 fsFlash.setSaturation(phaseIdx, FsToolbox::value(fluidState.saturation(phaseIdx)));
359 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
360 fsFlash.setMoleFraction(phaseIdx, compIdx, FsToolbox::value(fluidState.moleFraction(phaseIdx, compIdx)));
363 paramCache.updateAll(fsFlash);
364 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
365 if (!FluidSystem::phaseIsActive(phaseIdx))
368 Scalar rho = FluidSystem::template density<FlashFluidState, Scalar>(fsFlash, paramCache, phaseIdx);
369 fsFlash.setDensity(phaseIdx, rho);
373 ComponentVector globalMolarities(0.0);
374 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
375 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
376 if (!FluidSystem::phaseIsActive(phaseIdx))
379 globalMolarities[compIdx] +=
380 fsFlash.saturation(phaseIdx) * fsFlash.molarity(phaseIdx, compIdx);
388 NcpFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
400 using ConstEvaluation =
typename std::remove_reference<typename FluidState::Scalar>::type;
401 using FsEvaluation =
typename std::remove_const<ConstEvaluation>::type;
402 using FsToolbox = MathToolbox<FsEvaluation>;
404 bool gasPresent = FluidSystem::phaseIsActive(gasPhaseIdx)?(fluidState.saturation(gasPhaseIdx) > 0.0):
false;
405 bool oilPresent = FluidSystem::phaseIsActive(oilPhaseIdx)?(fluidState.saturation(oilPhaseIdx) > 0.0):
false;
406 bool waterPresent = FluidSystem::phaseIsActive(waterPhaseIdx)?(fluidState.saturation(waterPhaseIdx) > 0.0):
false;
407 const auto& saltSaturation = BlackOil::getSaltSaturation_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
408 bool precipitatedSaltPresent = enableSaltPrecipitation?(saltSaturation > 0.0):
false;
409 bool oneActivePhases = FluidSystem::numActivePhases() == 1;
416 if (gasPresent && FluidSystem::enableVaporizedOil() && !oilPresent){
417 primaryVarsMeaningPressure_ = PressureMeaning::Pg;
418 }
else if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
419 primaryVarsMeaningPressure_ = PressureMeaning::Po;
420 }
else if ( waterPresent && FluidSystem::enableDissolvedGasInWater() && !gasPresent){
421 primaryVarsMeaningPressure_ = PressureMeaning::Pw;
422 }
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
423 primaryVarsMeaningPressure_ = PressureMeaning::Pg;
425 assert(FluidSystem::phaseIsActive(waterPhaseIdx));
426 primaryVarsMeaningPressure_ = PressureMeaning::Pw;
433 if ( waterPresent && gasPresent ){
434 primaryVarsMeaningWater_ = WaterMeaning::Sw;
435 }
else if (gasPresent && FluidSystem::enableVaporizedWater()) {
436 primaryVarsMeaningWater_ = WaterMeaning::Rvw;
437 }
else if (waterPresent && FluidSystem::enableDissolvedGasInWater()) {
438 primaryVarsMeaningWater_ = WaterMeaning::Rsw;
439 }
else if (FluidSystem::phaseIsActive(waterPhaseIdx) && !oneActivePhases) {
440 primaryVarsMeaningWater_ = WaterMeaning::Sw;
442 primaryVarsMeaningWater_ = WaterMeaning::Disabled;
450 if ( gasPresent && oilPresent ) {
451 primaryVarsMeaningGas_ = GasMeaning::Sg;
452 }
else if (oilPresent && FluidSystem::enableDissolvedGas()) {
453 primaryVarsMeaningGas_ = GasMeaning::Rs;
454 }
else if (gasPresent && FluidSystem::enableVaporizedOil()){
455 primaryVarsMeaningGas_ = GasMeaning::Rv;
456 }
else if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx)) {
457 primaryVarsMeaningGas_ = GasMeaning::Sg;
459 primaryVarsMeaningGas_ = GasMeaning::Disabled;
463 if constexpr (enableSaltPrecipitation){
464 if (precipitatedSaltPresent)
465 primaryVarsMeaningBrine_ = BrineMeaning::Sp;
467 primaryVarsMeaningBrine_ = BrineMeaning::Cs;
469 primaryVarsMeaningBrine_ = BrineMeaning::Disabled;
474 case PressureMeaning::Po:
475 this->setScaledPressure_(FsToolbox::value(fluidState.pressure(oilPhaseIdx)));
477 case PressureMeaning::Pg:
478 this->setScaledPressure_(FsToolbox::value(fluidState.pressure(gasPhaseIdx)));
480 case PressureMeaning::Pw:
481 this->setScaledPressure_(FsToolbox::value(fluidState.pressure(waterPhaseIdx)));
484 throw std::logic_error(
"No valid primary variable selected for pressure");
487 case WaterMeaning::Sw:
489 (*this)[waterSwitchIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
492 case WaterMeaning::Rvw:
494 const auto& rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
495 (*this)[waterSwitchIdx] = rvw;
498 case WaterMeaning::Rsw:
500 const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
501 (*this)[waterSwitchIdx] = Rsw;
504 case WaterMeaning::Disabled:
509 throw std::logic_error(
"No valid primary variable selected for water");
514 (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
519 const auto& rs = BlackOil::getRs_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
520 (*this)[compositionSwitchIdx] = rs;
525 const auto& rv = BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
526 (*this)[compositionSwitchIdx] = rv;
529 case GasMeaning::Disabled:
534 throw std::logic_error(
"No valid primary variable selected for composision");
550 unsigned globalDofIdx,
551 [[maybe_unused]] Scalar swMaximum,
552 Scalar thresholdWaterFilledCell, Scalar eps = 0.0)
569 Scalar saltConcentration = 0.0;
570 const Scalar& T = asImp_().temperature_(problem, globalDofIdx);
572 sw = (*this)[waterSwitchIdx];
574 sg = (*this)[compositionSwitchIdx];
583 if constexpr (enableSaltPrecipitation) {
585 if (primaryVarsMeaningBrine() == BrineMeaning::Sp) {
586 saltConcentration = saltSolubility;
587 Scalar saltSat = (*this)[saltConcentrationIdx];
590 (*this)[saltConcentrationIdx] = saltSolubility;
593 else if (primaryVarsMeaningBrine() == BrineMeaning::Cs) {
594 saltConcentration = (*this)[saltConcentrationIdx];
595 if (saltConcentration > saltSolubility + eps){
597 (*this)[saltConcentrationIdx] = 0.0;
605 if constexpr (enableSolvent) {
606 if (SolventModule::isSolubleInWater()) {
607 Scalar p = (*this)[pressureSwitchIdx];
608 Scalar solLimit = SolventModule::solubilityLimit(
pvtRegionIndex(), T , p, saltConcentration);
609 if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
610 Scalar solSat = (*this)[solventSaturationIdx];
613 (*this)[solventSaturationIdx] = solLimit;
617 else if (primaryVarsMeaningSolvent() == SolventMeaning::Rsolw) {
618 Scalar rsolw = (*this)[solventSaturationIdx];
619 if (rsolw > solLimit + eps){
621 (*this)[solventSaturationIdx] = 0.0;
628 bool changed =
false;
636 if (sw >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {
639 if constexpr (waterEnabled) {
640 (*this)[Indices::waterSwitchIdx] = std::min(swMaximum, sw);
644 if constexpr (compositionSwitchEnabled)
645 (*
this)[Indices::compositionSwitchIdx] = 0.0;
649 if constexpr (compositionSwitchEnabled)
657 if (BrineModule::hasPcfactTables() && primaryVarsMeaningBrine() == BrineMeaning::Sp) {
658 unsigned satnumRegionIdx = problem.satnumRegionIndex(globalDofIdx);
659 Scalar Sp = saltConcentration_();
660 Scalar porosityFactor = min(1.0 - Sp, 1.0);
661 const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
662 pcFactor_ = pcfactTable.eval(porosityFactor,
true);
669 case WaterMeaning::Sw:
672 if(sw < -eps && sg > eps && FluidSystem::enableVaporizedWater()) {
673 Scalar p = this->pressure_();
675 std::array<Scalar, numPhases> pC = { 0.0 };
676 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
677 Scalar so = 1.0 - sg - solventSaturation_();
678 computeCapillaryPressures_(pC, so, sg + solventSaturation_(), 0.0, matParams);
679 p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
681 Scalar rvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
686 (*this)[Indices::waterSwitchIdx] = rvwSat;
692 if(sg < -eps && sw > eps && FluidSystem::enableDissolvedGasInWater()) {
693 const Scalar pg = this->pressure_();
695 std::array<Scalar, numPhases> pC = { 0.0 };
696 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
697 Scalar so = 1.0 - sw - solventSaturation_();
698 computeCapillaryPressures_(pC, so, 0.0, sw, matParams);
699 Scalar pw = pg + pcFactor_ * (pC[waterPhaseIdx] - pC[gasPhaseIdx]);
700 Scalar rswSat = FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
705 Scalar rswMax = problem.maxGasDissolutionFactor(0, globalDofIdx);
706 (*this)[Indices::waterSwitchIdx] = min(rswSat, rswMax);
708 this->setScaledPressure_(pw);
714 case WaterMeaning::Rvw:
716 const Scalar& rvw = (*this)[waterSwitchIdx];
717 Scalar p = this->pressure_();
719 std::array<Scalar, numPhases> pC = { 0.0 };
720 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
721 Scalar so = 1.0 - sg - solventSaturation_();
722 computeCapillaryPressures_(pC, so, sg + solventSaturation_(), 0.0, matParams);
723 p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
725 Scalar rvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
730 if (rvw > rvwSat*(1.0 + eps)) {
732 (*this)[Indices::waterSwitchIdx] = 0.0;
737 case WaterMeaning::Rsw:
742 const Scalar& pw = this->pressure_();
744 Scalar rswSat = FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
749 Scalar rsw = (*this)[Indices::waterSwitchIdx];
750 Scalar rswMax = problem.maxGasDissolutionFactor(0, globalDofIdx);
751 if (rsw > min(rswSat, rswMax)) {
754 (*this)[Indices::waterSwitchIdx] = 1.0;
756 std::array<Scalar, numPhases> pC = { 0.0 };
757 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
758 computeCapillaryPressures_(pC, 0.0, 0.0, 1.0, matParams);
759 Scalar pg = pw + pcFactor_ * (pC[gasPhaseIdx] - pC[waterPhaseIdx]);
760 this->setScaledPressure_(pg);
765 case WaterMeaning::Disabled:
770 throw std::logic_error(
"No valid primary variable selected for water");
784 Scalar s = 1.0 - sw - solventSaturation_();
785 if (sg < -eps && s > 0.0 && FluidSystem::enableDissolvedGas()) {
786 const Scalar po = this->pressure_();
788 Scalar soMax = std::max(s, problem.maxOilSaturation(globalDofIdx));
789 Scalar rsMax = problem.maxGasDissolutionFactor(0, globalDofIdx);
793 : FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
798 (*this)[Indices::compositionSwitchIdx] = std::min(rsMax, rsSat);
801 Scalar so = 1.0 - sw - solventSaturation_() - sg;
802 if (so < -eps && sg > 0.0 && FluidSystem::enableVaporizedOil()) {
807 const Scalar po = this->pressure_();
808 std::array<Scalar, numPhases> pC = { 0.0 };
809 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
810 computeCapillaryPressures_(pC, 0.0, sg + solventSaturation_(), sw, matParams);
811 Scalar pg = po + pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
816 this->setScaledPressure_(pg);
817 Scalar soMax = problem.maxOilSaturation(globalDofIdx);
818 Scalar rvMax = problem.maxOilVaporizationFactor(0, globalDofIdx);
822 : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
828 (*this)[Indices::compositionSwitchIdx] = std::min(rvMax, rvSat);
838 const Scalar po = this->pressure_();
839 Scalar so = 1.0 - sw - solventSaturation_();
840 Scalar soMax = std::max(so, problem.maxOilSaturation(globalDofIdx));
841 Scalar rsMax = problem.maxGasDissolutionFactor(0, globalDofIdx);
845 : FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
851 Scalar rs = (*this)[Indices::compositionSwitchIdx];
852 if (rs > std::min(rsMax, rsSat*(Scalar{1.0} + eps))) {
855 (*this)[Indices::compositionSwitchIdx] = 0.0;
866 const Scalar pg = this->pressure_();
867 Scalar soMax = problem.maxOilSaturation(globalDofIdx);
868 Scalar rvMax = problem.maxOilVaporizationFactor(0, globalDofIdx);
872 : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
878 Scalar rv = (*this)[Indices::compositionSwitchIdx];
879 if (rv > std::min(rvMax, rvSat*(Scalar{1.0} + eps))) {
883 Scalar sg2 = 1.0 - sw - solventSaturation_();
884 std::array<Scalar, numPhases> pC = { 0.0 };
885 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
886 computeCapillaryPressures_(pC,
888 sg2 + solventSaturation_(),
891 Scalar po = pg + pcFactor_ * (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
895 this->setScaledPressure_(po);
896 (*this)[Indices::compositionSwitchIdx] = sg2;
901 case GasMeaning::Disabled:
906 throw std::logic_error(
"No valid primary variable selected for water");