My Project
Loading...
Searching...
No Matches
blackoilintensivequantities.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_HH
29#define EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_HH
30
31#include "blackoilproperties.hh"
42
43#include <opm/common/TimingMacros.hpp>
44#include <opm/common/OpmLog/OpmLog.hpp>
45
46#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
47
48#include <opm/material/fluidstates/BlackOilFluidState.hpp>
49#include <opm/material/common/Valgrind.hpp>
50
52
53#include <opm/utility/CopyablePtr.hpp>
54
55#include <dune/common/fmatrix.hh>
56
57#include <cstring>
58#include <utility>
59
60namespace Opm {
68template <class TypeTag>
70 : public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
71 , public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
72 , public BlackOilDiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
73 , public BlackOilDispersionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDispersion>() >
75 , public BlackOilExtboIntensiveQuantities<TypeTag>
77 , public BlackOilFoamIntensiveQuantities<TypeTag>
78 , public BlackOilBrineIntensiveQuantities<TypeTag>
79 , public BlackOilEnergyIntensiveQuantities<TypeTag>
80 , public BlackOilMICPIntensiveQuantities<TypeTag>
82{
85
95
97 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
99 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
100 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
101 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
102 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
103 enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
104 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
105 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
106 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
107 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
108 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
109 enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
110 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
112 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
113 enum { waterCompIdx = FluidSystem::waterCompIdx };
114 enum { oilCompIdx = FluidSystem::oilCompIdx };
115 enum { gasCompIdx = FluidSystem::gasCompIdx };
116 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
117 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
118 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
119 enum { dimWorld = GridView::dimensionworld };
120 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
121
122 static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
123 static constexpr bool waterEnabled = Indices::waterEnabled;
124 static constexpr bool gasEnabled = Indices::gasEnabled;
125 static constexpr bool oilEnabled = Indices::oilEnabled;
126
127 using Toolbox = MathToolbox<Evaluation>;
128 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
129 using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
132
133 using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
137
138
139public:
140 using FluidState = BlackOilFluidState<Evaluation,
141 FluidSystem,
142 enableTemperature,
143 enableEnergy,
144 compositionSwitchEnabled,
145 enableVapwat,
146 enableBrine,
147 enableSaltPrecipitation,
148 has_disgas_in_water,
149 Indices::numPhases>;
150 using ScalarFluidState = BlackOilFluidState<Scalar,
151 FluidSystem,
152 enableTemperature,
153 enableEnergy,
154 compositionSwitchEnabled,
155 enableVapwat,
156 enableBrine,
157 enableSaltPrecipitation,
158 has_disgas_in_water,
159 Indices::numPhases>;
161
163 {
164 if (compositionSwitchEnabled) {
165 fluidState_.setRs(0.0);
166 fluidState_.setRv(0.0);
167 }
168 if (enableVapwat) {
169 fluidState_.setRvw(0.0);
170 }
171 if (has_disgas_in_water) {
172 fluidState_.setRsw(0.0);
173 }
174 }
176
177 BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default;
178
179 void updateTempSalt(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
180 {
181 if constexpr (enableTemperature || enableEnergy) {
182 asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
183 }
184
185 if constexpr (enableBrine) {
186 asImp_().updateSaltConcentration_(elemCtx, dofIdx, timeIdx);
187 }
188 }
189
190 void updateSaturations(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
191 {
192 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
193
194 // extract the water and the gas saturations for convenience
195 Evaluation Sw = 0.0;
196 if constexpr (waterEnabled) {
197 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
198 assert(Indices::waterSwitchIdx >= 0);
199 if constexpr (Indices::waterSwitchIdx >= 0) {
200 Sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
201 }
202 } else if(priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw ||
203 priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Disabled) {
204 // water is enabled but is not a primary variable i.e. one component/phase case
205 // or two-phase water + gas with only water present
206 Sw = 1.0;
207 } // else i.e. for MeaningWater() = Rvw, Sw is still 0.0;
208 }
209 Evaluation Sg = 0.0;
210 if constexpr (gasEnabled) {
211 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
212 assert(Indices::compositionSwitchIdx >= 0);
213 if constexpr (compositionSwitchEnabled) {
214 Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
215 }
216 } else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
217 Sg = 1.0 - Sw;
218 } else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Disabled) {
219 if constexpr (waterEnabled) {
220 Sg = 1.0 - Sw; // two phase water + gas
221 } else {
222 // one phase case
223 Sg = 1.0;
224 }
225 }
226 }
227 Valgrind::CheckDefined(Sg);
228 Valgrind::CheckDefined(Sw);
229
230 Evaluation So = 1.0 - Sw - Sg;
231
232 // deal with solvent
233 if constexpr (enableSolvent) {
234 if(priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
235 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
236 So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
237 } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
238 Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
239 }
240 }
241 }
242
243 if (FluidSystem::phaseIsActive(waterPhaseIdx))
244 fluidState_.setSaturation(waterPhaseIdx, Sw);
245
246 if (FluidSystem::phaseIsActive(gasPhaseIdx))
247 fluidState_.setSaturation(gasPhaseIdx, Sg);
248
249 if (FluidSystem::phaseIsActive(oilPhaseIdx))
250 fluidState_.setSaturation(oilPhaseIdx, So);
251 }
252
253 void updateRelpermAndPressures(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
254 {
255 const auto& problem = elemCtx.problem();
256 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
257 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
258
259 // Solvent saturation manipulation:
260 // After this, gas saturation will actually be (gas sat + solvent sat)
261 // until set back to just gas saturation in the corresponding call to
262 // solventPostSatFuncUpdate_() further down.
263 if constexpr (enableSolvent) {
264 asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
265 }
266
267 // Phase relperms.
268 problem.updateRelperms(mobility_, dirMob_, fluidState_, globalSpaceIdx);
269
270 // now we compute all phase pressures
271 std::array<Evaluation, numPhases> pC;
272 const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
273 MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
274
275 // scaling the capillary pressure due to salt precipitation
276 if constexpr (enableBrine) {
277 if (BrineModule::hasPcfactTables() && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
278 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, dofIdx, timeIdx);
279 const Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
280 const Evaluation porosityFactor = min(1.0 - Sp, 1.0); //phi/phi_0
281 const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
282 const Evaluation pcFactor = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
283 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
284 if (FluidSystem::phaseIsActive(phaseIdx)) {
285 pC[phaseIdx] *= pcFactor;
286 }
287 }
288 }
289
290 // oil is the reference phase for pressure
291 if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
292 const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
293 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
294 if (FluidSystem::phaseIsActive(phaseIdx))
295 fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
296 } else if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pw) {
297 const Evaluation& pw = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
298 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
299 if (FluidSystem::phaseIsActive(phaseIdx))
300 fluidState_.setPressure(phaseIdx, pw + (pC[phaseIdx] - pC[waterPhaseIdx]));
301 } else {
302 assert(FluidSystem::phaseIsActive(oilPhaseIdx));
303 const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
304 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
305 if (FluidSystem::phaseIsActive(phaseIdx))
306 fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
307 }
308
309 // Update the Saturation functions for the blackoil solvent module.
310 // Including setting gas saturation back to hydrocarbon gas saturation.
311 // Note that this depend on the pressures, so it must be called AFTER the pressures
312 // have been updated.
313 if constexpr (enableSolvent) {
314 asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
315 }
316
317 }
318
319 Evaluation updateRsRvRsw(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
320 {
321 const auto& problem = elemCtx.problem();
322 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
323 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
324 const unsigned pvtRegionIdx = priVars.pvtRegionIndex();
325
326 Scalar RvMax = FluidSystem::enableVaporizedOil()
327 ? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
328 : 0.0;
329 Scalar RsMax = FluidSystem::enableDissolvedGas()
330 ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
331 : 0.0;
332 Scalar RswMax = FluidSystem::enableDissolvedGasInWater()
333 ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
334 : 0.0;
335
336 Evaluation SoMax = 0.0;
337 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
338 SoMax = max(fluidState_.saturation(oilPhaseIdx),
339 problem.maxOilSaturation(globalSpaceIdx));
340 }
341 // take the meaning of the switching primary variable into account for the gas
342 // and oil phase compositions
343 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
344 const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
345 fluidState_.setRs(Rs);
346 } else {
347 if (FluidSystem::enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
348 const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
349 FluidSystem::saturatedDissolutionFactor(fluidState_,
350 oilPhaseIdx,
351 pvtRegionIdx,
352 SoMax);
353 fluidState_.setRs(min(RsMax, RsSat));
354 }
355 else if constexpr (compositionSwitchEnabled)
356 fluidState_.setRs(0.0);
357 }
358 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
359 const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
360 fluidState_.setRv(Rv);
361 } else {
362 if (FluidSystem::enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
363 const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
364 FluidSystem::saturatedDissolutionFactor(fluidState_,
365 gasPhaseIdx,
366 pvtRegionIdx,
367 SoMax);
368 fluidState_.setRv(min(RvMax, RvSat));
369 }
370 else if constexpr (compositionSwitchEnabled)
371 fluidState_.setRv(0.0);
372 }
373
374 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rvw) {
375 const auto& Rvw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
376 fluidState_.setRvw(Rvw);
377 } else {
378 if (FluidSystem::enableVaporizedWater()) { // Add Sg > 0? i.e. if only water set rv = 0)
379 const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
380 gasPhaseIdx,
381 pvtRegionIdx);
382 fluidState_.setRvw(RvwSat);
383 }
384 }
385
386 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw) {
387 const auto& Rsw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
388 fluidState_.setRsw(Rsw);
389 } else {
390 if (FluidSystem::enableDissolvedGasInWater()) {
391 const Evaluation& RswSat = FluidSystem::saturatedDissolutionFactor(fluidState_,
392 waterPhaseIdx,
393 pvtRegionIdx);
394 fluidState_.setRsw(min(RswMax, RswSat));
395 }
396 }
397
398 return SoMax;
399 }
400
401 void updateMobilityAndInvB()
402 {
403 const unsigned pvtRegionIdx = fluidState_.pvtRegionIndex();
404
405 // compute the phase densities and transform the phase permeabilities into mobilities
406 int nmobilities = 1;
407 std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
408 if (dirMob_) {
409 for (int i=0; i<3; i++) {
410 nmobilities += 1;
411 mobilities.push_back(&(dirMob_->getArray(i)));
412 }
413 }
414 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
415 if (!FluidSystem::phaseIsActive(phaseIdx))
416 continue;
417 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
418 fluidState_.setInvB(phaseIdx, b);
419 const auto& mu = FluidSystem::viscosity(fluidState_, phaseIdx, pvtRegionIdx);
420 for (int i = 0; i<nmobilities; i++) {
421 if (enableExtbo && phaseIdx == oilPhaseIdx) {
422 (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
423 }
424 else if (enableExtbo && phaseIdx == gasPhaseIdx) {
425 (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
426 }
427 else {
428 (*mobilities[i])[phaseIdx] /= mu;
429 }
430 }
431 }
432 Valgrind::CheckDefined(mobility_);
433 }
434
435 void updatePhaseDensities()
436 {
437 const unsigned pvtRegionIdx = fluidState_.pvtRegionIndex();
438
439 // calculate the phase densities
440 Evaluation rho;
441 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
442 rho = fluidState_.invB(waterPhaseIdx);
443 rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
444 if (FluidSystem::enableDissolvedGasInWater()) {
445 rho +=
446 fluidState_.invB(waterPhaseIdx) *
447 fluidState_.Rsw() *
448 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
449 }
450 fluidState_.setDensity(waterPhaseIdx, rho);
451 }
452
453 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
454 rho = fluidState_.invB(gasPhaseIdx);
455 rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
456 if (FluidSystem::enableVaporizedOil()) {
457 rho +=
458 fluidState_.invB(gasPhaseIdx) *
459 fluidState_.Rv() *
460 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
461 }
462 if (FluidSystem::enableVaporizedWater()) {
463 rho +=
464 fluidState_.invB(gasPhaseIdx) *
465 fluidState_.Rvw() *
466 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
467 }
468 fluidState_.setDensity(gasPhaseIdx, rho);
469 }
470
471 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
472 rho = fluidState_.invB(oilPhaseIdx);
473 rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
474 if (FluidSystem::enableDissolvedGas()) {
475 rho +=
476 fluidState_.invB(oilPhaseIdx) *
477 fluidState_.Rs() *
478 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
479 }
480 fluidState_.setDensity(oilPhaseIdx, rho);
481 }
482 }
483
484 void updatePorosity(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
485 {
486 const auto& problem = elemCtx.problem();
487 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
488 const auto& linearizationType = problem.model().linearizer().getLinearizationType();
489 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
490
491 // retrieve the porosity from the problem
492 referencePorosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
493 porosity_ = referencePorosity_;
494
495 // the porosity must be modified by the compressibility of the
496 // rock...
497 Scalar rockCompressibility = problem.rockCompressibility(globalSpaceIdx);
498 if (rockCompressibility > 0.0) {
499 Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
500 Evaluation x;
501 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
502 x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
503 } else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
504 x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
505 } else {
506 x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
507 }
508 porosity_ *= 1.0 + x + 0.5*x*x;
509 }
510
511 // deal with water induced rock compaction
512 porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
513
514 // deal with MICP
515 if constexpr (enableMICP){
516 Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
517 Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
518 // minimum porosity of 1e-8 to prevent numerical issues
519 porosity_ -= min(biofilm_ + calcite_, referencePorosity_ - 1e-8);
520 }
521
522 // deal with salt-precipitation
523 if (enableSaltPrecipitation && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
524 Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
525 porosity_ *= (1.0 - Sp);
526 }
527 }
528
529 void assertFiniteMembers()
530 {
531 // some safety checks in debug mode
532 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
533 if (!FluidSystem::phaseIsActive(phaseIdx))
534 continue;
535
536 assert(isfinite(fluidState_.density(phaseIdx)));
537 assert(isfinite(fluidState_.saturation(phaseIdx)));
538 assert(isfinite(fluidState_.temperature(phaseIdx)));
539 assert(isfinite(fluidState_.pressure(phaseIdx)));
540 assert(isfinite(fluidState_.invB(phaseIdx)));
541 }
542 assert(isfinite(fluidState_.Rs()));
543 assert(isfinite(fluidState_.Rv()));
544 }
545
549 void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
550 {
551 ParentType::update(elemCtx, dofIdx, timeIdx);
552
554
555 const auto& problem = elemCtx.problem();
556 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
557 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
558 const unsigned pvtRegionIdx = priVars.pvtRegionIndex();
559
560 fluidState_.setPvtRegionIndex(pvtRegionIdx);
561
562 updateTempSalt(elemCtx, dofIdx, timeIdx);
563 updateSaturations(elemCtx, dofIdx, timeIdx);
564 updateRelpermAndPressures(elemCtx, dofIdx, timeIdx);
565
566 // update extBO parameters
567 if constexpr (enableExtbo) {
568 asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
569 }
570
571 Evaluation SoMax = updateRsRvRsw(elemCtx, dofIdx, timeIdx);
572
573 updateMobilityAndInvB();
574 updatePhaseDensities();
575 updatePorosity(elemCtx, dofIdx, timeIdx);
576
577 rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*this, globalSpaceIdx);
578
579 if constexpr (enableSolvent) {
580 asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
581 }
582 if constexpr (enableExtbo) {
583 asImp_().zPvtUpdate_();
584 }
585 if constexpr (enablePolymer) {
586 asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
587 }
588
589 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
590 paramCache.setRegionIndex(pvtRegionIdx);
591 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
592 paramCache.setMaxOilSat(SoMax);
593 }
594 paramCache.updateAll(fluidState_);
595
596 if constexpr (enableEnergy) {
597 asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
598 }
599 if constexpr (enableFoam) {
600 asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
601 }
602 if constexpr (enableMICP) {
603 asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
604 }
605 if constexpr (enableBrine) {
606 asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
607 }
608 if constexpr (enableConvectiveMixing) {
609 // The ifs are here is to avoid extra calculations for
610 // cases without CO2STORE and DRSDTCON.
611 if (problem.simulator().vanguard().eclState().runspec().co2Storage()) {
612 if (problem.drsdtconIsActive(globalSpaceIdx, problem.simulator().episodeIndex())) {
613 asImp_().updateSaturatedDissolutionFactor_();
614 }
615 }
616 }
617
618 // update the quantities which are required by the chosen
619 // velocity model
620 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
621
622 // update the diffusion specific quantities of the intensive quantities
623 DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
624
625 // update the dispersion specific quantities of the intensive quantities
626 DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
627
628#ifndef NDEBUG
629 assertFiniteMembers();
630#endif
631 }
632
636 const FluidState& fluidState() const
637 { return fluidState_; }
638
642 const Evaluation& mobility(unsigned phaseIdx) const
643 { return mobility_[phaseIdx]; }
644
645 const Evaluation& mobility(unsigned phaseIdx, FaceDir::DirEnum facedir) const
646 {
647 using Dir = FaceDir::DirEnum;
648 if (dirMob_) {
649 switch(facedir) {
650 case Dir::XMinus:
651 case Dir::XPlus:
652 return dirMob_->mobilityX_[phaseIdx];
653 case Dir::YMinus:
654 case Dir::YPlus:
655 return dirMob_->mobilityY_[phaseIdx];
656 case Dir::ZMinus:
657 case Dir::ZPlus:
658 return dirMob_->mobilityZ_[phaseIdx];
659 default:
660 throw std::runtime_error("Unexpected face direction");
661 }
662 }
663 else {
664 return mobility_[phaseIdx];
665 }
666
667 }
668
672 const Evaluation& porosity() const
673 { return porosity_; }
674
678 const Evaluation& rockCompTransMultiplier() const
679 { return rockCompTransMultiplier_; }
680
689 -> decltype(std::declval<FluidState>().pvtRegionIndex())
690 { return fluidState_.pvtRegionIndex(); }
691
695 Evaluation relativePermeability(unsigned phaseIdx) const
696 {
697 // warning: slow
698 return fluidState_.viscosity(phaseIdx)*mobility(phaseIdx);
699 }
700
707 Scalar referencePorosity() const
708 { return referencePorosity_; }
709
710 const Evaluation& permFactor() const
711 {
712 if constexpr (enableMICP) {
713 return MICPIntQua::permFactor();
714 }
715 else if constexpr (enableSaltPrecipitation) {
716 return BrineIntQua::permFactor();
717 }
718 else {
719 throw std::logic_error("permFactor() called but salt precipitation or MICP are disabled");
720 }
721 }
722
723private:
724 friend BlackOilSolventIntensiveQuantities<TypeTag>;
725 friend BlackOilExtboIntensiveQuantities<TypeTag>;
726 friend BlackOilPolymerIntensiveQuantities<TypeTag>;
727 friend BlackOilEnergyIntensiveQuantities<TypeTag>;
728 friend BlackOilFoamIntensiveQuantities<TypeTag>;
730 friend BlackOilMICPIntensiveQuantities<TypeTag>;
731
732 Implementation& asImp_()
733 { return *static_cast<Implementation*>(this); }
734
735 FluidState fluidState_;
736 Scalar referencePorosity_;
737 Evaluation porosity_;
738 Evaluation rockCompTransMultiplier_;
739 std::array<Evaluation,numPhases> mobility_;
740
741 // Instead of writing a custom copy constructor and a custom assignment operator just to handle
742 // the dirMob_ unique ptr member variable when copying BlackOilIntensiveQuantites (see for example
743 // updateIntensitiveQuantities_() in fvbaseelementcontext.hh for a copy example) we write the below
744 // custom wrapper class CopyablePtr which wraps the unique ptr and makes it copyable.
745 //
746 // The advantage of this approach is that we avoid having to call all the base class copy constructors and
747 // assignment operators explicitly (which is needed when writing the custom copy constructor and assignment
748 // operators) which could become a maintenance burden. For example, when adding a new base class (if that should
749 // be needed sometime in the future) to BlackOilIntensiveQuantites we could forget to update the copy
750 // constructor and assignment operators.
751 //
752 // We want each copy of the BlackOilIntensiveQuantites to be unique, (TODO: why?) so we have to make a copy
753 // of the unique_ptr each time we copy construct or assign to it from another BlackOilIntensiveQuantites.
754 // (On the other hand, if a copy could share the ptr with the original, a shared_ptr could be used instead and the
755 // wrapper would not be needed)
756 DirectionalMobilityPtr dirMob_;
757};
758
759} // namespace Opm
760
761#endif
Contains the classes required to extend the black-oil model by brine.
Classes required for dynamic convective mixing.
Classes required for molecular diffusion.
Classes required for mechanical dispersion.
Contains the classes required to extend the black-oil model by energy.
Contains the classes required to extend the black-oil model by solvent component.
Contains the classes required to extend the black-oil model to include the effects of foam.
Contains the classes required to extend the black-oil model by MICP.
Contains the classes required to extend the black-oil model by polymer.
Declares the properties required by the black oil model.
Contains the classes required to extend the black-oil model by solvents.
Provides the volumetric quantities required for the equations needed by the brine extension of the bl...
Definition blackoilbrinemodules.hh:340
Contains the high level supplements required to extend the black oil model by brine.
Definition blackoilbrinemodules.hh:50
Provides the volumetric quantities required for the equations needed by the convective mixing (DRSDTC...
Definition blackoilconvectivemixingmodule.hh:387
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes.
Definition blackoildiffusionmodule.hh:309
Provides the volumetric quantities required for the calculation of dispersive fluxes.
Definition blackoildispersionmodule.hh:308
Provides the volumetric quantities required for the equations needed by the energys extension of the ...
Definition blackoilenergymodules.hh:333
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition blackoilextbomodules.hh:399
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition blackoilfoammodules.hh:381
Contains the quantities which are are constant within a finite volume in the black-oil model.
Definition blackoilintensivequantities.hh:82
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition blackoilintensivequantities.hh:672
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition blackoilintensivequantities.hh:642
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition blackoilintensivequantities.hh:549
Evaluation relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition blackoilintensivequantities.hh:695
auto pvtRegionIndex() const -> decltype(std::declval< FluidState >().pvtRegionIndex())
Returns the index of the PVT region used to calculate the thermodynamic quantities.
Definition blackoilintensivequantities.hh:688
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition blackoilintensivequantities.hh:636
const Evaluation & rockCompTransMultiplier() const
The pressure-dependent transmissibility multiplier due to rock compressibility.
Definition blackoilintensivequantities.hh:678
Scalar referencePorosity() const
Returns the porosity of the rock at reference conditions.
Definition blackoilintensivequantities.hh:707
Provides the volumetric quantities required for the equations needed by the MICP extension of the bla...
Definition blackoilmicpmodules.hh:410
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition blackoilpolymermodules.hh:552
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition blackoilsolventmodules.hh:520
This file contains definitions related to directional mobilities.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
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