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