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"
41
42#include <opm/common/TimingMacros.hpp>
43#include <opm/common/OpmLog/OpmLog.hpp>
44
45#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
46
47#include <opm/material/fluidstates/BlackOilFluidState.hpp>
48#include <opm/material/common/Valgrind.hpp>
49
51
52#include <opm/utility/CopyablePtr.hpp>
53
54#include <dune/common/fmatrix.hh>
55
56#include <cstring>
57#include <utility>
58
59namespace Opm {
67template <class TypeTag>
69 : public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
70 , public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
71 , public BlackOilDiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
72 , public BlackOilDispersionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDispersion>() >
74 , public BlackOilExtboIntensiveQuantities<TypeTag>
76 , public BlackOilFoamIntensiveQuantities<TypeTag>
77 , public BlackOilBrineIntensiveQuantities<TypeTag>
78 , public BlackOilEnergyIntensiveQuantities<TypeTag>
79 , public BlackOilMICPIntensiveQuantities<TypeTag>
80{
83
93
95 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
97 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
100 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
101 enum { has_disgas_in_water = getPropValue<TypeTag, Properties::EnableDisgasInWater>() };
102 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
103 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
104 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
105 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
106 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
107 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
109 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
110 enum { waterCompIdx = FluidSystem::waterCompIdx };
111 enum { oilCompIdx = FluidSystem::oilCompIdx };
112 enum { gasCompIdx = FluidSystem::gasCompIdx };
113 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
114 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
115 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
116 enum { dimWorld = GridView::dimensionworld };
117 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
118
119 static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
120 static constexpr bool waterEnabled = Indices::waterEnabled;
121 static constexpr bool gasEnabled = Indices::gasEnabled;
122 static constexpr bool oilEnabled = Indices::oilEnabled;
123
124 using Toolbox = MathToolbox<Evaluation>;
125 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
126 using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
129
130 using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
132
133
134public:
135 using FluidState = BlackOilFluidState<Evaluation,
136 FluidSystem,
137 enableTemperature,
138 enableEnergy,
139 compositionSwitchEnabled,
140 enableVapwat,
141 enableBrine,
142 enableSaltPrecipitation,
143 has_disgas_in_water,
144 Indices::numPhases>;
145 using ScalarFluidState = BlackOilFluidState<Scalar,
146 FluidSystem,
147 enableTemperature,
148 enableEnergy,
149 compositionSwitchEnabled,
150 enableVapwat,
151 enableBrine,
152 enableSaltPrecipitation,
153 has_disgas_in_water,
154 Indices::numPhases>;
156
158 {
159 if (compositionSwitchEnabled) {
160 fluidState_.setRs(0.0);
161 fluidState_.setRv(0.0);
162 }
163 if (enableVapwat) {
164 fluidState_.setRvw(0.0);
165 }
166 if (has_disgas_in_water) {
167 fluidState_.setRsw(0.0);
168 }
169 }
171
172 BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default;
173
174 void updateTempSalt(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
175 {
176 if constexpr (enableTemperature || enableEnergy) {
177 asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
178 }
179
180 if constexpr (enableBrine) {
181 asImp_().updateSaltConcentration_(elemCtx, dofIdx, timeIdx);
182 }
183 }
184
185 void updateSaturations(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
186 {
187 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
188
189 // extract the water and the gas saturations for convenience
190 Evaluation Sw = 0.0;
191 if constexpr (waterEnabled) {
192 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
193 assert(Indices::waterSwitchIdx >= 0);
194 if constexpr (Indices::waterSwitchIdx >= 0) {
195 Sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
196 }
197 } else if(priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw ||
198 priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Disabled) {
199 // water is enabled but is not a primary variable i.e. one component/phase case
200 // or two-phase water + gas with only water present
201 Sw = 1.0;
202 } // else i.e. for MeaningWater() = Rvw, Sw is still 0.0;
203 }
204 Evaluation Sg = 0.0;
205 if constexpr (gasEnabled) {
206 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
207 assert(Indices::compositionSwitchIdx >= 0);
208 if constexpr (compositionSwitchEnabled) {
209 Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
210 }
211 } else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
212 Sg = 1.0 - Sw;
213 } else if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Disabled) {
214 if constexpr (waterEnabled) {
215 Sg = 1.0 - Sw; // two phase water + gas
216 } else {
217 // one phase case
218 Sg = 1.0;
219 }
220 }
221 }
222 Valgrind::CheckDefined(Sg);
223 Valgrind::CheckDefined(Sw);
224
225 Evaluation So = 1.0 - Sw - Sg;
226
227 // deal with solvent
228 if constexpr (enableSolvent) {
229 if(priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
230 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
231 So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
232 } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
233 Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
234 }
235 }
236 }
237
238 if (FluidSystem::phaseIsActive(waterPhaseIdx))
239 fluidState_.setSaturation(waterPhaseIdx, Sw);
240
241 if (FluidSystem::phaseIsActive(gasPhaseIdx))
242 fluidState_.setSaturation(gasPhaseIdx, Sg);
243
244 if (FluidSystem::phaseIsActive(oilPhaseIdx))
245 fluidState_.setSaturation(oilPhaseIdx, So);
246 }
247
248 void updateRelpermAndPressures(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
249 {
250 const auto& problem = elemCtx.problem();
251 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
252 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
253
254 // Solvent saturation manipulation:
255 // After this, gas saturation will actually be (gas sat + solvent sat)
256 // until set back to just gas saturation in the corresponding call to
257 // solventPostSatFuncUpdate_() further down.
258 if constexpr (enableSolvent) {
259 asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
260 }
261
262 // Phase relperms.
263 problem.updateRelperms(mobility_, dirMob_, fluidState_, globalSpaceIdx);
264
265 // now we compute all phase pressures
266 std::array<Evaluation, numPhases> pC;
267 const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
268 MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
269
270 // scaling the capillary pressure due to salt precipitation
271 if constexpr (enableBrine) {
272 if (BrineModule::hasPcfactTables() && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
273 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, dofIdx, timeIdx);
274 const Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
275 const Evaluation porosityFactor = min(1.0 - Sp, 1.0); //phi/phi_0
276 const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
277 const Evaluation pcFactor = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
278 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
279 if (FluidSystem::phaseIsActive(phaseIdx)) {
280 pC[phaseIdx] *= pcFactor;
281 }
282 }
283 }
284
285 // oil is the reference phase for pressure
286 if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pg) {
287 const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
288 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
289 if (FluidSystem::phaseIsActive(phaseIdx))
290 fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
291 } else if (priVars.primaryVarsMeaningPressure() == PrimaryVariables::PressureMeaning::Pw) {
292 const Evaluation& pw = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
293 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
294 if (FluidSystem::phaseIsActive(phaseIdx))
295 fluidState_.setPressure(phaseIdx, pw + (pC[phaseIdx] - pC[waterPhaseIdx]));
296 } else {
297 assert(FluidSystem::phaseIsActive(oilPhaseIdx));
298 const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
299 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
300 if (FluidSystem::phaseIsActive(phaseIdx))
301 fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
302 }
303
304 // Update the Saturation functions for the blackoil solvent module.
305 // Including setting gas saturation back to hydrocarbon gas saturation.
306 // Note that this depend on the pressures, so it must be called AFTER the pressures
307 // have been updated.
308 if constexpr (enableSolvent) {
309 asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
310 }
311
312 }
313
314 Evaluation updateRsRvRsw(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
315 {
316 const auto& problem = elemCtx.problem();
317 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
318 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
319 const unsigned pvtRegionIdx = priVars.pvtRegionIndex();
320
321 Scalar RvMax = FluidSystem::enableVaporizedOil()
322 ? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
323 : 0.0;
324 Scalar RsMax = FluidSystem::enableDissolvedGas()
325 ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
326 : 0.0;
327 Scalar RswMax = FluidSystem::enableDissolvedGasInWater()
328 ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
329 : 0.0;
330
331 Evaluation SoMax = 0.0;
332 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
333 SoMax = max(fluidState_.saturation(oilPhaseIdx),
334 problem.maxOilSaturation(globalSpaceIdx));
335 }
336 // take the meaning of the switching primary variable into account for the gas
337 // and oil phase compositions
338 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
339 const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
340 fluidState_.setRs(Rs);
341 } else {
342 if (FluidSystem::enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
343 const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
344 FluidSystem::saturatedDissolutionFactor(fluidState_,
345 oilPhaseIdx,
346 pvtRegionIdx,
347 SoMax);
348 fluidState_.setRs(min(RsMax, RsSat));
349 }
350 else if constexpr (compositionSwitchEnabled)
351 fluidState_.setRs(0.0);
352 }
353 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
354 const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
355 fluidState_.setRv(Rv);
356 } else {
357 if (FluidSystem::enableVaporizedOil() ) { // Add Sg > 0? i.e. if only water set rv = 0)
358 const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
359 FluidSystem::saturatedDissolutionFactor(fluidState_,
360 gasPhaseIdx,
361 pvtRegionIdx,
362 SoMax);
363 fluidState_.setRv(min(RvMax, RvSat));
364 }
365 else if constexpr (compositionSwitchEnabled)
366 fluidState_.setRv(0.0);
367 }
368
369 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rvw) {
370 const auto& Rvw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
371 fluidState_.setRvw(Rvw);
372 } else {
373 if (FluidSystem::enableVaporizedWater()) { // Add Sg > 0? i.e. if only water set rv = 0)
374 const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
375 gasPhaseIdx,
376 pvtRegionIdx);
377 fluidState_.setRvw(RvwSat);
378 }
379 }
380
381 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Rsw) {
382 const auto& Rsw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx);
383 fluidState_.setRsw(Rsw);
384 } else {
385 if (FluidSystem::enableDissolvedGasInWater()) {
386 const Evaluation& RswSat = FluidSystem::saturatedDissolutionFactor(fluidState_,
387 waterPhaseIdx,
388 pvtRegionIdx);
389 fluidState_.setRsw(min(RswMax, RswSat));
390 }
391 }
392
393 return SoMax;
394 }
395
396 void updateMobilityAndInvB()
397 {
398 const unsigned pvtRegionIdx = fluidState_.pvtRegionIndex();
399
400 // compute the phase densities and transform the phase permeabilities into mobilities
401 int nmobilities = 1;
402 std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
403 if (dirMob_) {
404 for (int i=0; i<3; i++) {
405 nmobilities += 1;
406 mobilities.push_back(&(dirMob_->getArray(i)));
407 }
408 }
409 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
410 if (!FluidSystem::phaseIsActive(phaseIdx))
411 continue;
412 const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
413 fluidState_.setInvB(phaseIdx, b);
414 const auto& mu = FluidSystem::viscosity(fluidState_, phaseIdx, pvtRegionIdx);
415 for (int i = 0; i<nmobilities; i++) {
416 if (enableExtbo && phaseIdx == oilPhaseIdx) {
417 (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
418 }
419 else if (enableExtbo && phaseIdx == gasPhaseIdx) {
420 (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
421 }
422 else {
423 (*mobilities[i])[phaseIdx] /= mu;
424 }
425 }
426 }
427 Valgrind::CheckDefined(mobility_);
428 }
429
430 void updatePhaseDensities()
431 {
432 const unsigned pvtRegionIdx = fluidState_.pvtRegionIndex();
433
434 // calculate the phase densities
435 Evaluation rho;
436 if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
437 rho = fluidState_.invB(waterPhaseIdx);
438 rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
439 if (FluidSystem::enableDissolvedGasInWater()) {
440 rho +=
441 fluidState_.invB(waterPhaseIdx) *
442 fluidState_.Rsw() *
443 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
444 }
445 fluidState_.setDensity(waterPhaseIdx, rho);
446 }
447
448 if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
449 rho = fluidState_.invB(gasPhaseIdx);
450 rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
451 if (FluidSystem::enableVaporizedOil()) {
452 rho +=
453 fluidState_.invB(gasPhaseIdx) *
454 fluidState_.Rv() *
455 FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
456 }
457 if (FluidSystem::enableVaporizedWater()) {
458 rho +=
459 fluidState_.invB(gasPhaseIdx) *
460 fluidState_.Rvw() *
461 FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
462 }
463 fluidState_.setDensity(gasPhaseIdx, rho);
464 }
465
466 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
467 rho = fluidState_.invB(oilPhaseIdx);
468 rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
469 if (FluidSystem::enableDissolvedGas()) {
470 rho +=
471 fluidState_.invB(oilPhaseIdx) *
472 fluidState_.Rs() *
473 FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
474 }
475 fluidState_.setDensity(oilPhaseIdx, rho);
476 }
477 }
478
479 void updatePorosity(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
480 {
481 const auto& problem = elemCtx.problem();
482 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
483 const auto& linearizationType = problem.model().linearizer().getLinearizationType();
484 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
485
486 // retrieve the porosity from the problem
487 referencePorosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
488 porosity_ = referencePorosity_;
489
490 // the porosity must be modified by the compressibility of the
491 // rock...
492 Scalar rockCompressibility = problem.rockCompressibility(globalSpaceIdx);
493 if (rockCompressibility > 0.0) {
494 Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
495 Evaluation x;
496 if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
497 x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
498 } else if (FluidSystem::phaseIsActive(waterPhaseIdx)){
499 x = rockCompressibility*(fluidState_.pressure(waterPhaseIdx) - rockRefPressure);
500 } else {
501 x = rockCompressibility*(fluidState_.pressure(gasPhaseIdx) - rockRefPressure);
502 }
503 porosity_ *= 1.0 + x + 0.5*x*x;
504 }
505
506 // deal with water induced rock compaction
507 porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
508
509 // the MICP processes change the porosity
510 if constexpr (enableMICP){
511 Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
512 Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
513 porosity_ += - biofilm_ - calcite_;
514 }
515
516 // deal with salt-precipitation
517 if (enableSaltPrecipitation && priVars.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
518 Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
519 porosity_ *= (1.0 - Sp);
520 }
521 }
522
523 void assertFiniteMembers()
524 {
525 // some safety checks in debug mode
526 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
527 if (!FluidSystem::phaseIsActive(phaseIdx))
528 continue;
529
530 assert(isfinite(fluidState_.density(phaseIdx)));
531 assert(isfinite(fluidState_.saturation(phaseIdx)));
532 assert(isfinite(fluidState_.temperature(phaseIdx)));
533 assert(isfinite(fluidState_.pressure(phaseIdx)));
534 assert(isfinite(fluidState_.invB(phaseIdx)));
535 }
536 assert(isfinite(fluidState_.Rs()));
537 assert(isfinite(fluidState_.Rv()));
538 }
539
543 void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
544 {
545 ParentType::update(elemCtx, dofIdx, timeIdx);
546
548
549 const auto& problem = elemCtx.problem();
550 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
551 const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
552 const unsigned pvtRegionIdx = priVars.pvtRegionIndex();
553
554 fluidState_.setPvtRegionIndex(pvtRegionIdx);
555
556 updateTempSalt(elemCtx, dofIdx, timeIdx);
557 updateSaturations(elemCtx, dofIdx, timeIdx);
558 updateRelpermAndPressures(elemCtx, dofIdx, timeIdx);
559
560 // update extBO parameters
561 if constexpr (enableExtbo) {
562 asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
563 }
564
565 Evaluation SoMax = updateRsRvRsw(elemCtx, dofIdx, timeIdx);
566
567 updateMobilityAndInvB();
568 updatePhaseDensities();
569 updatePorosity(elemCtx, dofIdx, timeIdx);
570
571 rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*this, globalSpaceIdx);
572
573 if constexpr (enableSolvent) {
574 asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
575 }
576 if constexpr (enableExtbo) {
577 asImp_().zPvtUpdate_();
578 }
579 if constexpr (enablePolymer) {
580 asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
581 }
582
583 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
584 paramCache.setRegionIndex(pvtRegionIdx);
585 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
586 paramCache.setMaxOilSat(SoMax);
587 }
588 paramCache.updateAll(fluidState_);
589
590 if constexpr (enableEnergy) {
591 asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
592 }
593 if constexpr (enableFoam) {
594 asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
595 }
596 if constexpr (enableMICP) {
597 asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
598 }
599 if constexpr (enableBrine) {
600 asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
601 }
602
603 // update the quantities which are required by the chosen
604 // velocity model
605 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
606
607 // update the diffusion specific quantities of the intensive quantities
608 DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
609
610 // update the dispersion specific quantities of the intensive quantities
611 DispersionIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
612
613#ifndef NDEBUG
614 assertFiniteMembers();
615#endif
616 }
617
621 const FluidState& fluidState() const
622 { return fluidState_; }
623
627 const Evaluation& mobility(unsigned phaseIdx) const
628 { return mobility_[phaseIdx]; }
629
630 const Evaluation& mobility(unsigned phaseIdx, FaceDir::DirEnum facedir) const
631 {
632 using Dir = FaceDir::DirEnum;
633 if (dirMob_) {
634 switch(facedir) {
635 case Dir::XMinus:
636 case Dir::XPlus:
637 return dirMob_->mobilityX_[phaseIdx];
638 case Dir::YMinus:
639 case Dir::YPlus:
640 return dirMob_->mobilityY_[phaseIdx];
641 case Dir::ZMinus:
642 case Dir::ZPlus:
643 return dirMob_->mobilityZ_[phaseIdx];
644 default:
645 throw std::runtime_error("Unexpected face direction");
646 }
647 }
648 else {
649 return mobility_[phaseIdx];
650 }
651
652 }
653
657 const Evaluation& porosity() const
658 { return porosity_; }
659
663 const Evaluation& rockCompTransMultiplier() const
664 { return rockCompTransMultiplier_; }
665
674 -> decltype(std::declval<FluidState>().pvtRegionIndex())
675 { return fluidState_.pvtRegionIndex(); }
676
680 Evaluation relativePermeability(unsigned phaseIdx) const
681 {
682 // warning: slow
683 return fluidState_.viscosity(phaseIdx)*mobility(phaseIdx);
684 }
685
692 Scalar referencePorosity() const
693 { return referencePorosity_; }
694
695private:
703
704 Implementation& asImp_()
705 { return *static_cast<Implementation*>(this); }
706
707 FluidState fluidState_;
708 Scalar referencePorosity_;
709 Evaluation porosity_;
710 Evaluation rockCompTransMultiplier_;
711 std::array<Evaluation,numPhases> mobility_;
712
713 // Instead of writing a custom copy constructor and a custom assignment operator just to handle
714 // the dirMob_ unique ptr member variable when copying BlackOilIntensiveQuantites (see for example
715 // updateIntensitiveQuantities_() in fvbaseelementcontext.hh for a copy example) we write the below
716 // custom wrapper class CopyablePtr which wraps the unique ptr and makes it copyable.
717 //
718 // The advantage of this approach is that we avoid having to call all the base class copy constructors and
719 // assignment operators explicitly (which is needed when writing the custom copy constructor and assignment
720 // operators) which could become a maintenance burden. For example, when adding a new base class (if that should
721 // be needed sometime in the future) to BlackOilIntensiveQuantites we could forget to update the copy
722 // constructor and assignment operators.
723 //
724 // We want each copy of the BlackOilIntensiveQuantites to be unique, (TODO: why?) so we have to make a copy
725 // of the unique_ptr each time we copy construct or assign to it from another BlackOilIntensiveQuantites.
726 // (On the other hand, if a copy could share the ptr with the original, a shared_ptr could be used instead and the
727 // wrapper would not be needed)
728 DirectionalMobilityPtr dirMob_;
729};
730
731} // namespace Opm
732
733#endif
Contains the classes required to extend the black-oil model by brine.
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.
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 calculation of molecular diffusive fluxes.
Definition blackoildiffusionmodule.hh:274
Provides the volumetric quantities required for the calculation of dispersive fluxes.
Definition blackoildispersionmodule.hh:283
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:80
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition blackoilintensivequantities.hh:657
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition blackoilintensivequantities.hh:627
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition blackoilintensivequantities.hh:543
Evaluation relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition blackoilintensivequantities.hh:680
auto pvtRegionIndex() const -> decltype(std::declval< FluidState >().pvtRegionIndex())
Returns the index of the PVT region used to calculate the thermodynamic quantities.
Definition blackoilintensivequantities.hh:673
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition blackoilintensivequantities.hh:621
const Evaluation & rockCompTransMultiplier() const
The pressure-dependent transmissibility multiplier due to rock compressibility.
Definition blackoilintensivequantities.hh:663
Scalar referencePorosity() const
Returns the porosity of the rock at reference conditions.
Definition blackoilintensivequantities.hh:692
Provides the volumetric quantities required for the equations needed by the MICP extension of the bla...
Definition blackoilmicpmodules.hh:371
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