179 void update(
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
181 ParentType::update(elemCtx, dofIdx, timeIdx);
182 OPM_TIMEBLOCK_LOCAL(blackoilIntensiveQuanititiesUpdate);
183 const auto& problem = elemCtx.problem();
184 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
185 const auto& linearizationType = problem.model().linearizer().getLinearizationType();
186 unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
187 Scalar RvMax = FluidSystem::enableVaporizedOil()
188 ? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
190 Scalar RsMax = FluidSystem::enableDissolvedGas()
191 ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
193 Scalar RswMax = FluidSystem::enableDissolvedGasInWater()
194 ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
197 asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
199 unsigned pvtRegionIdx = priVars.pvtRegionIndex();
200 fluidState_.setPvtRegionIndex(pvtRegionIdx);
202 asImp_().updateSaltConcentration_(elemCtx, dofIdx, timeIdx);
206 if constexpr (waterEnabled) {
207 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
208 Sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
209 }
else if(priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw ||
210 priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Disabled) {
217 if constexpr (gasEnabled) {
218 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
219 Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
220 }
else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
222 }
else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Disabled) {
223 if constexpr (waterEnabled) {
231 Valgrind::CheckDefined(Sg);
232 Valgrind::CheckDefined(Sw);
234 Evaluation So = 1.0 - Sw - Sg;
237 if constexpr (enableSolvent) {
238 if(priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
239 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
240 So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
241 }
else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
242 Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
247 if (FluidSystem::phaseIsActive(waterPhaseIdx))
248 fluidState_.setSaturation(waterPhaseIdx, Sw);
250 if (FluidSystem::phaseIsActive(gasPhaseIdx))
251 fluidState_.setSaturation(gasPhaseIdx, Sg);
253 if (FluidSystem::phaseIsActive(oilPhaseIdx))
254 fluidState_.setSaturation(oilPhaseIdx, So);
256 asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
259 std::array<Evaluation, numPhases> pC;
260 const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
261 MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
262 problem.updateRelperms(mobility_, dirMob_, fluidState_, globalSpaceIdx);
265 if (BrineModule::hasPcfactTables() && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
266 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, dofIdx, timeIdx);
267 const Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
268 const Evaluation porosityFactor = min(1.0 - Sp, 1.0);
269 const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
270 const Evaluation pcFactor = pcfactTable.eval(porosityFactor,
true);
271 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
272 if (FluidSystem::phaseIsActive(phaseIdx)) {
273 pC[phaseIdx] *= pcFactor;
278 if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
279 const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
280 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
281 if (FluidSystem::phaseIsActive(phaseIdx))
282 fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
283 }
else if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pw) {
284 const Evaluation& pw = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
285 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
286 if (FluidSystem::phaseIsActive(phaseIdx))
287 fluidState_.setPressure(phaseIdx, pw + (pC[phaseIdx] - pC[waterPhaseIdx]));
289 assert(FluidSystem::phaseIsActive(oilPhaseIdx));
290 const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
291 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
292 if (FluidSystem::phaseIsActive(phaseIdx))
293 fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
297 asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
300 asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
302 Evaluation SoMax = 0.0;
303 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
304 SoMax = max(fluidState_.saturation(oilPhaseIdx),
305 problem.maxOilSaturation(globalSpaceIdx));
309 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
310 const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
311 fluidState_.setRs(Rs);
313 if (FluidSystem::enableDissolvedGas()) {
314 const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
315 FluidSystem::saturatedDissolutionFactor(fluidState_,
319 fluidState_.setRs(min(RsMax, RsSat));
321 else if constexpr (compositionSwitchEnabled)
322 fluidState_.setRs(0.0);
324 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
325 const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
326 fluidState_.setRv(Rv);
328 if (FluidSystem::enableVaporizedOil() ) {
329 const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
330 FluidSystem::saturatedDissolutionFactor(fluidState_,
334 fluidState_.setRv(min(RvMax, RvSat));
336 else if constexpr (compositionSwitchEnabled)
337 fluidState_.setRv(0.0);
340 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rvw) {
341 const auto& Rvw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
342 fluidState_.setRvw(Rvw);
344 if (FluidSystem::enableVaporizedWater()) {
345 const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
348 fluidState_.setRvw(RvwSat);
352 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw) {
353 const auto& Rsw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
354 fluidState_.setRsw(Rsw);
356 if (FluidSystem::enableDissolvedGasInWater()) {
357 const Evaluation& RswSat = FluidSystem::saturatedDissolutionFactor(fluidState_,
360 fluidState_.setRsw(min(RswMax, RswSat));
364 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
365 paramCache.setRegionIndex(pvtRegionIdx);
366 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
367 paramCache.setMaxOilSat(SoMax);
369 paramCache.updateAll(fluidState_);
373 std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
375 for (
int i=0; i<3; i++) {
377 mobilities.push_back(&(dirMob_->getArray(i)));
380 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
381 if (!FluidSystem::phaseIsActive(phaseIdx))
383 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
384 fluidState_.setInvB(phaseIdx, b);
385 const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
386 for (
int i = 0; i<nmobilities; i++) {
387 if (enableExtbo && phaseIdx == oilPhaseIdx) {
388 (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
390 else if (enableExtbo && phaseIdx == gasPhaseIdx) {
391 (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
394 (*mobilities[i])[phaseIdx] /= mu;
398 Valgrind::CheckDefined(mobility_);
402 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
403 rho = fluidState_.invB(waterPhaseIdx);
404 rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
405 if (FluidSystem::enableDissolvedGasInWater()) {
407 fluidState_.invB(waterPhaseIdx) *
409 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
411 fluidState_.setDensity(waterPhaseIdx, rho);
414 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
415 rho = fluidState_.invB(gasPhaseIdx);
416 rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
417 if (FluidSystem::enableVaporizedOil()) {
419 fluidState_.invB(gasPhaseIdx) *
421 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
423 if (FluidSystem::enableVaporizedWater()) {
425 fluidState_.invB(gasPhaseIdx) *
427 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
429 fluidState_.setDensity(gasPhaseIdx, rho);
432 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
433 rho = fluidState_.invB(oilPhaseIdx);
434 rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
435 if (FluidSystem::enableDissolvedGas()) {
437 fluidState_.invB(oilPhaseIdx) *
439 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
441 fluidState_.setDensity(oilPhaseIdx, rho);
445 referencePorosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
446 porosity_ = referencePorosity_;
450 Scalar rockCompressibility = problem.rockCompressibility(globalSpaceIdx);
451 if (rockCompressibility > 0.0) {
452 Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
454 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
455 x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
456 }
else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
457 x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
459 x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
461 porosity_ *= 1.0 + x + 0.5*x*x;
465 porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*
this, globalSpaceIdx);
468 if constexpr (enableMICP){
469 Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
470 Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
471 porosity_ += - biofilm_ - calcite_;
475 if (enableSaltPrecipitation && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
476 Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
477 porosity_ *= (1.0 - Sp);
482 asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
483 asImp_().zPvtUpdate_();
484 asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
485 asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
486 asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
487 asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
488 asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
492 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
495 DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
498 DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
502 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
503 if (!FluidSystem::phaseIsActive(phaseIdx))
506 assert(isfinite(fluidState_.density(phaseIdx)));
507 assert(isfinite(fluidState_.saturation(phaseIdx)));
508 assert(isfinite(fluidState_.temperature(phaseIdx)));
509 assert(isfinite(fluidState_.pressure(phaseIdx)));
510 assert(isfinite(fluidState_.invB(phaseIdx)));
512 assert(isfinite(fluidState_.Rs()));
513 assert(isfinite(fluidState_.Rv()));