blackoilprimaryvariables.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 OPM is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 2 of the License, or
8 (at your option) any later version.
9 OPM is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with OPM. If not, see <http://www.gnu.org/licenses/>.
15 Consult the COPYING file in the top-level source directory of this
16 module for the precise wording of the license and the list of
17 copyright holders.
18*/
24#ifndef EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH
25#define EWOMS_BLACK_OIL_PRIMARY_VARIABLES_HH
26
27#include <dune/common/fvector.hh>
28
29#include <opm/material/common/MathToolbox.hpp>
30#include <opm/material/common/Valgrind.hpp>
31#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
32#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
33
39
41
42#include <algorithm>
43#include <array>
44#include <cassert>
45#include <cstddef>
46#include <cstdint>
47#include <stdexcept>
48#include <type_traits>
49
50namespace Opm::Parameters {
51
52template<class Scalar>
54{ static constexpr Scalar value = 1.0; };
55
56} // namespace Opm::Parameters
57
58namespace Opm {
59
65template <class TypeTag>
67{
70
78
79 // number of equations
80 enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
81
82 // primary variable indices
83 enum { waterSwitchIdx = Indices::waterSwitchIdx };
84 enum { pressureSwitchIdx = Indices::pressureSwitchIdx };
85 enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
86 enum { saltConcentrationIdx = Indices::saltConcentrationIdx };
87 enum { solventSaturationIdx = Indices::solventSaturationIdx };
88
89 static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
90 static constexpr bool waterEnabled = Indices::waterEnabled;
91 static constexpr bool gasEnabled = Indices::gasEnabled;
92 static constexpr bool oilEnabled = Indices::oilEnabled;
93
94 // phase indices from the fluid system
95 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
96 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
97 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
98 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
99
100 // component indices from the fluid system
101 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
102 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
103 enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
104 enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
105 enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
106 enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
107 enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
108 enum { enableVapwat = getPropValue<TypeTag, Properties::EnableVapwat>() };
109 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
110 enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
111 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
112 enum { gasCompIdx = FluidSystem::gasCompIdx };
113 enum { waterCompIdx = FluidSystem::waterCompIdx };
114 enum { oilCompIdx = FluidSystem::oilCompIdx };
115
116 using Toolbox = MathToolbox<Evaluation>;
117 using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
122
123 static_assert(numPhases == 3, "The black-oil model assumes three phases!");
124 static_assert(numComponents == 3, "The black-oil model assumes three components!");
125
126public:
127 enum class WaterMeaning : std::uint8_t {
128 Sw, // water saturation
129 Rvw, // vaporized water
130 Rsw, // dissolved gas in water
131 Disabled, // The primary variable is not used
132 };
133
134 enum class PressureMeaning : std::uint8_t {
135 Po, // oil pressure
136 Pg, // gas pressure
137 Pw, // water pressure
138 };
139
140 enum class GasMeaning : std::uint8_t {
141 Sg, // gas saturation
142 Rs, // dissolved gas in oil
143 Rv, // vapporized oil
144 Disabled, // The primary variable is not used
145 };
146
147 enum class BrineMeaning : std::uint8_t {
148 Cs, // salt concentration
149 Sp, // (precipitated) salt saturation
150 Disabled, // The primary variable is not used
151 };
152
153 enum class SolventMeaning : std::uint8_t {
154 Ss, // solvent saturation
155 Rsolw, // dissolved solvent in water
156 Disabled, // The primary variable is not used
157 };
158
160 {
161 Valgrind::SetUndefined(*this);
162 pvtRegionIdx_ = 0;
163 }
164
169
171 {
173 result.pvtRegionIdx_ = 1;
174 result.primaryVarsMeaningBrine_ = BrineMeaning::Sp;
175 result.primaryVarsMeaningGas_ = GasMeaning::Rv;
176 result.primaryVarsMeaningPressure_ = PressureMeaning::Pg;
177 result.primaryVarsMeaningWater_ = WaterMeaning::Rsw;
178 result.primaryVarsMeaningSolvent_ = SolventMeaning::Ss;
179 for (std::size_t i = 0; i < result.size(); ++i) {
180 result[i] = i + 1;
181 }
182
183 return result;
184 }
185
186 static void init()
187 {
188 // TODO: these parameters have undocumented non-trivial dependencies
189 pressureScale_ = Parameters::Get<Parameters::PressureScale<Scalar>>();
190 }
191
192 static void registerParameters()
193 {
194 Parameters::Register<Parameters::PressureScale<Scalar>>
195 ("Scaling of pressure primary variable");
196 }
197
198 void setPressureScale(Scalar val)
199 { pressureScale_ = val; }
200
201 Evaluation
202 makeEvaluation(unsigned varIdx, unsigned timeIdx,
203 LinearizationType linearizationType = LinearizationType()) const
204 {
205 const Scalar scale = varIdx == pressureSwitchIdx ? this->pressureScale_ : Scalar{1.0};
206 if (std::is_same_v<Evaluation, Scalar>) {
207 return (*this)[varIdx] * scale; // finite differences
208 }
209 else {
210 // automatic differentiation
211 if (timeIdx == linearizationType.time) {
212 return Toolbox::createVariable((*this)[varIdx], varIdx) * scale;
213 }
214 else {
215 return Toolbox::createConstant((*this)[varIdx]) * scale;
216 }
217 }
218 }
219
227 void setPvtRegionIndex(unsigned value)
228 { pvtRegionIdx_ = static_cast<unsigned short>(value); }
229
233 unsigned pvtRegionIndex() const
234 { return pvtRegionIdx_; }
235
241 { return primaryVarsMeaningWater_; }
242
248 { primaryVarsMeaningWater_ = newMeaning; }
249
255 { return primaryVarsMeaningPressure_; }
256
262 { primaryVarsMeaningPressure_ = newMeaning; }
263
269 { return primaryVarsMeaningGas_; }
270
276 { primaryVarsMeaningGas_ = newMeaning; }
277
279 { return primaryVarsMeaningBrine_; }
280
286 { primaryVarsMeaningBrine_ = newMeaning; }
287
289 { return primaryVarsMeaningSolvent_; }
290
296 { primaryVarsMeaningSolvent_ = newMeaning; }
297
301 template <class FluidState>
302 void assignNaive(const FluidState& fluidState)
303 {
304 using ConstEvaluation = std::remove_reference_t<typename FluidState::Scalar>;
305 using FsEvaluation = std::remove_const_t<ConstEvaluation>;
306 using FsToolbox = MathToolbox<FsEvaluation>;
307
308 const bool gasPresent =
309 FluidSystem::phaseIsActive(gasPhaseIdx)
310 ? fluidState.saturation(gasPhaseIdx) > 0.0
311 : false;
312 const bool oilPresent =
313 FluidSystem::phaseIsActive(oilPhaseIdx)
314 ? fluidState.saturation(oilPhaseIdx) > 0.0
315 : false;
316 const bool waterPresent =
317 FluidSystem::phaseIsActive(waterPhaseIdx)
318 ? fluidState.saturation(waterPhaseIdx) > 0.0
319 : false;
320 const auto& saltSaturation =
321 BlackOil::getSaltSaturation_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
322 const bool precipitatedSaltPresent = enableSaltPrecipitation ? saltSaturation > 0.0 : false;
323 const bool oneActivePhases = FluidSystem::numActivePhases() == 1;
324 // deal with the primary variables for the energy extension
325 EnergyModule::assignPrimaryVars(*this, fluidState);
326
327 // Determine the meaning of the pressure primary variables
328 // Depending on the phases present, this variable is either interpreted as the
329 // pressure of the oil phase, gas phase (if no oil) or water phase (if only water)
330 if (gasPresent && FluidSystem::enableVaporizedOil() && !oilPresent) {
331 primaryVarsMeaningPressure_ = PressureMeaning::Pg;
332 }
333 else if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
334 primaryVarsMeaningPressure_ = PressureMeaning::Po;
335 }
336 else if (waterPresent && FluidSystem::enableDissolvedGasInWater() && !gasPresent) {
337 primaryVarsMeaningPressure_ = PressureMeaning::Pw;
338 }
339 else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
340 primaryVarsMeaningPressure_ = PressureMeaning::Pg;
341 }
342 else {
343 assert(FluidSystem::phaseIsActive(waterPhaseIdx));
344 primaryVarsMeaningPressure_ = PressureMeaning::Pw;
345 }
346
347 // Determine the meaning of the water primary variables
348 // Depending on the phases present, this variable is either interpreted as
349 // water saturation or vapporized water in the gas phase
350 // For two-phase gas-oil models and one-phase case the variable is disabled.
351 if (waterPresent && gasPresent) {
352 primaryVarsMeaningWater_ = WaterMeaning::Sw;
353 }
354 else if (gasPresent && FluidSystem::enableVaporizedWater()) {
355 primaryVarsMeaningWater_ = WaterMeaning::Rvw;
356 }
357 else if (waterPresent && FluidSystem::enableDissolvedGasInWater()) {
358 primaryVarsMeaningWater_ = WaterMeaning::Rsw;
359 }
360 else if (FluidSystem::phaseIsActive(waterPhaseIdx) && !oneActivePhases) {
361 primaryVarsMeaningWater_ = WaterMeaning::Sw;
362 }
363 else {
364 primaryVarsMeaningWater_ = WaterMeaning::Disabled;
365 }
366
367 // Determine the meaning of the gas primary variables
368 // Depending on the phases present, this variable is either interpreted as the
369 // saturation of the gas phase, as the fraction of the gas component in the oil
370 // phase (Rs) or as the fraction of the oil component (Rv) in the gas phase.
371 // For two-phase water-oil and water-gas models and one-phase case the variable is disabled.
372 if (gasPresent && oilPresent) {
373 primaryVarsMeaningGas_ = GasMeaning::Sg;
374 }
375 else if (oilPresent && FluidSystem::enableDissolvedGas()) {
376 primaryVarsMeaningGas_ = GasMeaning::Rs;
377 }
378 else if (gasPresent && FluidSystem::enableVaporizedOil()) {
379 primaryVarsMeaningGas_ = GasMeaning::Rv;
380 }
381 else if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx)) {
382 primaryVarsMeaningGas_ = GasMeaning::Sg;
383 }
384 else {
385 primaryVarsMeaningGas_ = GasMeaning::Disabled;
386 }
387
388 // Determine the meaning of the brine primary variables
389 if constexpr (enableSaltPrecipitation) {
390 if (precipitatedSaltPresent) {
391 primaryVarsMeaningBrine_ = BrineMeaning::Sp;
392 }
393 else {
394 primaryVarsMeaningBrine_ = BrineMeaning::Cs;
395 }
396 }
397 else {
398 primaryVarsMeaningBrine_ = BrineMeaning::Disabled;
399 }
400
401 // assign the actual primary variables
402 switch (primaryVarsMeaningPressure()) {
404 this->setScaledPressure_(FsToolbox::value(fluidState.pressure(oilPhaseIdx)));
405 break;
407 this->setScaledPressure_(FsToolbox::value(fluidState.pressure(gasPhaseIdx)));
408 break;
410 this->setScaledPressure_(FsToolbox::value(fluidState.pressure(waterPhaseIdx)));
411 break;
412 default:
413 throw std::logic_error("No valid primary variable selected for pressure");
414 }
415
416 switch (primaryVarsMeaningWater()) {
417 case WaterMeaning::Sw:
418 {
419 (*this)[waterSwitchIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
420 break;
421 }
423 {
424 const auto& rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
425 (*this)[waterSwitchIdx] = rvw;
426 break;
427 }
429 {
430 const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
431 (*this)[waterSwitchIdx] = Rsw;
432 break;
433 }
435 break;
436 default:
437 throw std::logic_error("No valid primary variable selected for water");
438 }
439 switch (primaryVarsMeaningGas()) {
440 case GasMeaning::Sg:
441 (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
442 break;
443 case GasMeaning::Rs:
444 {
445 const auto& rs = BlackOil::getRs_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
446 (*this)[compositionSwitchIdx] = rs;
447 break;
448 }
449 case GasMeaning::Rv:
450 {
451 const auto& rv = BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
452 (*this)[compositionSwitchIdx] = rv;
453 break;
454 }
456 break;
457 default:
458 throw std::logic_error("No valid primary variable selected for composision");
459 }
460 }
461
473 bool adaptPrimaryVariables(const Problem& problem,
474 unsigned globalDofIdx,
475 [[maybe_unused]] Scalar swMaximum,
476 Scalar thresholdWaterFilledCell, Scalar eps = 0.0)
477 {
478 // this function accesses quite a few black-oil specific low-level functions
479 // directly for better performance (instead of going the canonical way through
480 // the IntensiveQuantities). The reason is that most intensive quantities are not
481 // required to be able to decide if the primary variables needs to be switched or
482 // not, so it would be a waste to compute them.
483
484 // Both the primary variable meaning of water and gas are disabled i.e.
485 // It is a one-phase case and we no variable meaning switch is needed.
488 {
489 return false;
490 }
491
492 // Read the current saturation from the primary variables
493 Scalar sw = 0.0;
494 Scalar sg = 0.0;
495 Scalar saltConcentration = 0.0;
496 const Scalar& T = asImp_().temperature_(problem, globalDofIdx);
498 sw = (*this)[waterSwitchIdx];
499 }
501 sg = (*this)[compositionSwitchIdx];
502 }
504 sw = 1.0;
505 }
506
507 if (primaryVarsMeaningGas() == GasMeaning::Disabled && gasEnabled) {
508 sg = 1.0 - sw; // water + gas case
509 }
510
511 // if solid phase disappeares: Sp (Solid salt saturation) -> Cs (salt concentration)
512 // if solid phase appears: Cs (salt concentration) -> Sp (Solid salt saturation)
513 if constexpr (enableSaltPrecipitation) {
514 const Scalar saltSolubility = BrineModule::saltSol(pvtRegionIndex());
516 saltConcentration = saltSolubility;
517 const Scalar saltSat = (*this)[saltConcentrationIdx];
518 if (saltSat < -eps) { // precipitated salt dissappears
520 (*this)[saltConcentrationIdx] = saltSolubility; // set salt concentration to solubility limit
521 }
522 }
524 saltConcentration = (*this)[saltConcentrationIdx];
525 if (saltConcentration > saltSolubility + eps) { // salt concentration exceeds solubility limit
527 (*this)[saltConcentrationIdx] = 0.0;
528 }
529 }
530 }
531
532 // if solvent saturation disappeares: Ss (Solvent saturation) -> Rsolw (solvent dissolved in water)
533 // if solvent saturation appears: Rsolw (solvent dissolved in water) -> Ss (Solvent saturation)
534 // Scalar rsolw = 0.0; // not needed at the moment since we dont allow for vapwat in combination with rsolw
535 if constexpr (enableSolvent) {
537 const Scalar p = (*this)[pressureSwitchIdx]; // cap-pressure?
538 const Scalar solLimit =
539 SolventModule::solubilityLimit(pvtRegionIndex(), T , p, saltConcentration);
541 const Scalar solSat = (*this)[solventSaturationIdx];
542 if (solSat < -eps) { // solvent dissappears
544 (*this)[solventSaturationIdx] = solLimit; // set rsolw to solubility limit
545 }
546 }
548 const Scalar rsolw = (*this)[solventSaturationIdx];
549 if (rsolw > solLimit + eps) { // solvent appears as phase
551 (*this)[solventSaturationIdx] = 0.0;
552 }
553 }
554 }
555 }
556
557 // keep track if any meaning has changed
558 bool changed = false;
559
560 // Special case for cells with almost only water
561 // for these cells both saturations (if the phase is enabled) is used
562 // to avoid singular systems.
563 // If dissolved gas in water is enabled we shouldn't enter
564 // here but instead switch to Rsw as primary variable
565 // as sw >= 1.0 -> gas <= 0 (i.e. gas phase disappears)
566 if (sw >= thresholdWaterFilledCell && !FluidSystem::enableDissolvedGasInWater()) {
567 // make sure water saturations does not exceed sw_maximum. Default to 1.0
568 if constexpr (waterEnabled) {
569 (*this)[Indices::waterSwitchIdx] = std::min(swMaximum, sw);
571 }
572 // the hydrocarbon gas saturation is set to 0.0
573 if constexpr (compositionSwitchEnabled) {
574 (*this)[Indices::compositionSwitchIdx] = 0.0;
575 }
576
578 if (changed) {
579 if constexpr (compositionSwitchEnabled) {
581 }
582
583 // use water pressure?
584 }
585 return changed;
586 }
587
589 const unsigned satnumRegionIdx = problem.satnumRegionIndex(globalDofIdx);
590 const Scalar Sp = saltConcentration_();
591 const Scalar porosityFactor = std::min(1.0 - Sp, 1.0); //phi/phi_0
592 const auto& pcfactTable = BrineModule::pcfactTable(satnumRegionIdx);
593 pcFactor_ = pcfactTable.eval(porosityFactor, /*extrapolation=*/true);
594 }
595 else {
596 pcFactor_ = 1.0;
597 }
598
599 switch (primaryVarsMeaningWater()) {
600 case WaterMeaning::Sw:
601 {
602 // if water phase disappeares: Sw (water saturation) -> Rvw (fraction of water in gas phase)
603 if (sw < -eps && sg > eps && FluidSystem::enableVaporizedWater()) {
604 Scalar p = this->pressure_();
606 std::array<Scalar, numPhases> pC{};
607 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
608 const Scalar so = 1.0 - sg - solventSaturation_();
609 computeCapillaryPressures_(pC, so, sg + solventSaturation_(), /*sw=*/ 0.0, matParams);
610 p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
611 }
612 const Scalar rvwSat =
613 FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
614 T,
615 p,
616 saltConcentration);
618 (*this)[Indices::waterSwitchIdx] = rvwSat; // primary variable becomes Rvw
619 changed = true;
620 break;
621 }
622 // if gas phase disappeares: Sw (water saturation) -> Rsw (fraction of gas in water phase)
623 // and Pg (gas pressure) -> Pw ( water pressure)
624 if (sg < -eps && sw > eps && FluidSystem::enableDissolvedGasInWater()) {
625 const Scalar pg = this->pressure_();
627 std::array<Scalar, numPhases> pC = { 0.0 };
628 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
629 const Scalar so = 1.0 - sw - solventSaturation_();
630 computeCapillaryPressures_(pC, so, /*sg=*/ 0.0, sw, matParams);
631 const Scalar pw = pg + pcFactor_ * (pC[waterPhaseIdx] - pC[gasPhaseIdx]);
632 const Scalar rswSat =
633 FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
634 T,
635 pw,
636 saltConcentration);
638 const Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
639 (*this)[Indices::waterSwitchIdx] = std::min(rswSat, rswMax); //primary variable becomes Rsw
641 this->setScaledPressure_(pw);
642 changed = true;
643 break;
644 }
645 break;
646 }
648 {
649 const Scalar& rvw = (*this)[waterSwitchIdx];
650 Scalar p = this->pressure_();
652 std::array<Scalar, numPhases> pC{};
653 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
654 const Scalar so = 1.0 - sg - solventSaturation_();
655 computeCapillaryPressures_(pC, so, sg + solventSaturation_(), /*sw=*/ 0.0, matParams);
656 p += pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
657 }
658 const Scalar rvwSat =
659 FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
660 T,
661 p,
662 saltConcentration);
663 // if water phase appears: Rvw (fraction of water in gas phase) -> Sw (water saturation)
664 if (rvw > rvwSat * (1.0 + eps)) {
666 (*this)[Indices::waterSwitchIdx] = 0.0; // water saturation
667 changed = true;
668 }
669 break;
670 }
672 {
673 // Gas phase not present. The hydrocarbon gas phase
674 // appears as soon as more of the gas component is present in the water phase
675 // than what saturated water can hold.
676 const Scalar& pw = this->pressure_();
678 const Scalar rswSat =
679 FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
680 T,
681 pw,
682 saltConcentration);
683
684 const Scalar rsw = (*this)[Indices::waterSwitchIdx];
685 const Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
686 if (rsw > std::min(rswSat, rswMax)) {
687 // the gas phase appears, i.e., switch the primary variables to WaterMeaning::Sw
689 (*this)[Indices::waterSwitchIdx] = 1.0; // hydrocarbon water saturation
691 std::array<Scalar, numPhases> pC{};
692 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
693 computeCapillaryPressures_(pC, /*so=*/ 0.0, /*sg=*/ 0.0, /*sw=*/ 1.0, matParams);
694 const Scalar pg = pw + pcFactor_ * (pC[gasPhaseIdx] - pC[waterPhaseIdx]);
695 this->setScaledPressure_(pg);
696 changed = true;
697 }
698 break;
699 }
701 break;
702 default:
703 throw std::logic_error("No valid primary variable selected for water");
704 }
705
706 // if gas phase disappeares: Sg (gas saturation) -> Rs (fraction of gas in oil phase)
707 // if oil phase disappeares: Sg (gas saturation) -> Rv (fraction of oil in gas phase)
708 // Po (oil pressure ) -> Pg (gas pressure)
709
710 // if gas phase appears: Rs (fraction of gas in oil phase) -> Sg (gas saturation)
711 // if oil phase appears: Rv (fraction of oil in gas phase) -> Sg (gas saturation)
712 // Pg (gas pressure ) -> Po (oil pressure)
713 switch (primaryVarsMeaningGas()) {
714 case GasMeaning::Sg:
715 {
716 const Scalar s = 1.0 - sw - solventSaturation_();
717 if (sg < -eps && s > 0.0 && FluidSystem::enableDissolvedGas()) {
718 const Scalar po = this->pressure_();
720 const Scalar soMax = std::max(s, problem.maxOilSaturation(globalDofIdx));
721 const Scalar rsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
722 const Scalar rsSat =
723 enableExtbo
724 ? ExtboModule::rs(pvtRegionIndex(), po, zFraction_())
725 : FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
726 T,
727 po,
728 s,
729 soMax);
730 (*this)[Indices::compositionSwitchIdx] = std::min(rsMax, rsSat);
731 changed = true;
732 }
733 const Scalar so = 1.0 - sw - solventSaturation_() - sg;
734 if (so < -eps && sg > 0.0 && FluidSystem::enableVaporizedOil()) {
735 // the oil phase disappeared and some hydrocarbon gas phase is still
736 // present, i.e., switch the primary variables to GasMeaning::Rv.
737 // we only have the oil pressure readily available, but we need the gas
738 // pressure, i.e. we must determine capillary pressure
739 const Scalar po = this->pressure_();
740 std::array<Scalar, numPhases> pC{};
741 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
742 computeCapillaryPressures_(pC, /*so=*/0.0, sg + solventSaturation_(), sw, matParams);
743 const Scalar pg = po + pcFactor_ * (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
744
745 // we start at the GasMeaning::Rv value that corresponds to that of oil-saturated
746 // hydrocarbon gas
748 this->setScaledPressure_(pg);
749 const Scalar soMax = problem.maxOilSaturation(globalDofIdx);
750 const Scalar rvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
751 const Scalar rvSat =
752 enableExtbo
753 ? ExtboModule::rv(pvtRegionIndex(), pg, zFraction_())
754 : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
755 T,
756 pg,
757 Scalar(0),
758 soMax);
760 (*this)[Indices::compositionSwitchIdx] = std::min(rvMax, rvSat);
761 changed = true;
762 }
763 break;
764 }
765 case GasMeaning::Rs:
766 {
767 // Gas phase not present. The hydrocarbon gas phase
768 // appears as soon as more of the gas component is present in the oil phase
769 // than what saturated oil can hold.
770 const Scalar po = this->pressure_();
771 const Scalar so = 1.0 - sw - solventSaturation_();
772 const Scalar soMax = std::max(so, problem.maxOilSaturation(globalDofIdx));
773 const Scalar rsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
774 const Scalar rsSat =
775 enableExtbo
776 ? ExtboModule::rs(pvtRegionIndex(), po, zFraction_())
777 : FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
778 T,
779 po,
780 so,
781 soMax);
782
783 const Scalar rs = (*this)[Indices::compositionSwitchIdx];
784 if (rs > std::min(rsMax, rsSat * (Scalar{1.0} + eps))) {
785 // the gas phase appears, i.e., switch the primary variables to GasMeaning::Sg
787 (*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
788 changed = true;
789 }
790 break;
791 }
792 case GasMeaning::Rv:
793 {
794 // The oil phase appears as
795 // soon as more of the oil component is present in the hydrocarbon gas phase
796 // than what saturated gas contains. Note that we use the blackoil specific
797 // low-level PVT objects here for performance reasons.
798 const Scalar pg = this->pressure_();
799 const Scalar soMax = problem.maxOilSaturation(globalDofIdx);
800 const Scalar rvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
801 const Scalar rvSat =
802 enableExtbo
803 ? ExtboModule::rv(pvtRegionIndex(), pg, zFraction_())
804 : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
805 T,
806 pg,
807 /*so=*/Scalar(0.0),
808 soMax);
809
810 const Scalar rv = (*this)[Indices::compositionSwitchIdx];
811 if (rv > std::min(rvMax, rvSat * (Scalar{1.0} + eps))) {
812 // switch to phase equilibrium mode because the oil phase appears. here
813 // we also need the capillary pressures to calculate the oil phase
814 // pressure using the gas phase pressure
815 const Scalar sg2 = 1.0 - sw - solventSaturation_();
816 std::array<Scalar, numPhases> pC{};
817 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
818 computeCapillaryPressures_(pC,
819 /*so=*/0.0,
820 /*sg=*/sg2 + solventSaturation_(),
821 sw,
822 matParams);
823 const Scalar po = pg + pcFactor_ * (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
824
827 this->setScaledPressure_(po);
828 (*this)[Indices::compositionSwitchIdx] = sg2; // hydrocarbon gas saturation
829 changed = true;
830 }
831 break;
832 }
834 break;
835 default:
836 throw std::logic_error("No valid primary variable selected for water");
837 }
838 return changed;
839 }
840
842 {
845 {
846 return false;
847 }
848 Scalar sw = 0.0;
850 sw = (*this)[Indices::waterSwitchIdx];
851 }
852 Scalar sg = 0.0;
854 sg = (*this)[Indices::compositionSwitchIdx];
855 }
856
857 Scalar ssol = 0.0;
859 ssol =(*this) [Indices::solventSaturationIdx];
860 }
861
862 Scalar so = 1.0 - sw - sg - ssol;
863 sw = std::min(std::max(sw, Scalar{0.0}), Scalar{1.0});
864 so = std::min(std::max(so, Scalar{0.0}), Scalar{1.0});
865 sg = std::min(std::max(sg, Scalar{0.0}), Scalar{1.0});
866 ssol = std::min(std::max(ssol, Scalar{0.0}), Scalar{1.0});
867 const Scalar st = sw + so + sg + ssol;
868 sw = sw / st;
869 sg = sg / st;
870 ssol = ssol / st;
871 assert(st > 0.5);
873 (*this)[Indices::waterSwitchIdx] = sw;
874 }
876 (*this)[Indices::compositionSwitchIdx] = sg;
877 }
879 (*this) [Indices::solventSaturationIdx] = ssol;
880 }
881
882 return (st != 1);
883 }
884
886
887 using ParentType::operator=;
888
896 void checkDefined() const
897 {
898#ifndef NDEBUG
899 // check the "real" primary variables
900 for (unsigned i = 0; i < this->size(); ++i) {
901 Valgrind::CheckDefined((*this)[i]);
902 }
903
904 // check the "pseudo" primary variables
905 Valgrind::CheckDefined(primaryVarsMeaningWater_);
906 Valgrind::CheckDefined(primaryVarsMeaningGas_);
907 Valgrind::CheckDefined(primaryVarsMeaningPressure_);
908 Valgrind::CheckDefined(primaryVarsMeaningBrine_);
909 Valgrind::CheckDefined(primaryVarsMeaningSolvent_);
910
911 Valgrind::CheckDefined(pvtRegionIdx_);
912#endif // NDEBUG
913 }
914
915 template<class Serializer>
916 void serializeOp(Serializer& serializer)
917 {
918 using FV = Dune::FieldVector<Scalar, getPropValue<TypeTag, Properties::NumEq>()>;
919 serializer(static_cast<FV&>(*this));
920 serializer(primaryVarsMeaningWater_);
921 serializer(primaryVarsMeaningPressure_);
922 serializer(primaryVarsMeaningGas_);
923 serializer(primaryVarsMeaningBrine_);
924 serializer(primaryVarsMeaningSolvent_);
925 serializer(pvtRegionIdx_);
926 }
927
929 {
930 return
931 static_cast<const FvBasePrimaryVariables<TypeTag>&>(*this) == rhs
932 && this->primaryVarsMeaningWater_ == rhs.primaryVarsMeaningWater_
933 && this->primaryVarsMeaningPressure_ == rhs.primaryVarsMeaningPressure_
934 && this->primaryVarsMeaningGas_ == rhs.primaryVarsMeaningGas_
935 && this->primaryVarsMeaningBrine_ == rhs.primaryVarsMeaningBrine_
936 && this->primaryVarsMeaningSolvent_ == rhs.primaryVarsMeaningSolvent_
937 && this->pvtRegionIdx_ == rhs.pvtRegionIdx_;
938 }
939
940private:
941 Implementation& asImp_()
942 { return *static_cast<Implementation*>(this); }
943
944 const Implementation& asImp_() const
945 { return *static_cast<const Implementation*>(this); }
946
947 Scalar solventSaturation_() const
948 {
949 if constexpr (enableSolvent) {
951 return (*this)[Indices::solventSaturationIdx];
952 }
953 }
954 return 0.0;
955 }
956
957 Scalar zFraction_() const
958 {
959 if constexpr (enableExtbo) {
960 return (*this)[Indices::zFractionIdx];
961 }
962 else {
963 return 0.0;
964 }
965 }
966
967 Scalar saltConcentration_() const
968 {
969 if constexpr (enableBrine) {
970 return (*this)[Indices::saltConcentrationIdx];
971 }
972 else {
973 return 0.0;
974 }
975 }
976
977 Scalar temperature_(const Problem& problem, [[maybe_unused]] unsigned globalDofIdx) const
978 {
979 if constexpr (enableEnergy) {
980 return (*this)[Indices::temperatureIdx];
981 }
982 else if constexpr (enableTemperature) {
983 return problem.temperature(globalDofIdx, /*timeIdx*/ 0);
984 }
985 else {
986 return FluidSystem::reservoirTemperature();
987 }
988 }
989
990 template <class Container>
991 void computeCapillaryPressures_(Container& result,
992 Scalar so,
993 Scalar sg,
994 Scalar sw,
995 const MaterialLawParams& matParams) const
996 {
997 using SatOnlyFluidState = SimpleModularFluidState<Scalar,
998 numPhases,
999 numComponents,
1000 FluidSystem,
1001 /*storePressure=*/false,
1002 /*storeTemperature=*/false,
1003 /*storeComposition=*/false,
1004 /*storeFugacity=*/false,
1005 /*storeSaturation=*/true,
1006 /*storeDensity=*/false,
1007 /*storeViscosity=*/false,
1008 /*storeEnthalpy=*/false>;
1009 SatOnlyFluidState fluidState;
1010 fluidState.setSaturation(waterPhaseIdx, sw);
1011 fluidState.setSaturation(oilPhaseIdx, so);
1012 fluidState.setSaturation(gasPhaseIdx, sg);
1013
1014 MaterialLaw::capillaryPressures(result, matParams, fluidState);
1015 }
1016
1017 Scalar pressure_() const
1018 { return (*this)[Indices::pressureSwitchIdx] * this->pressureScale_; }
1019
1020 void setScaledPressure_(Scalar pressure)
1021 { (*this)[Indices::pressureSwitchIdx] = pressure / (this->pressureScale_); }
1022
1023 WaterMeaning primaryVarsMeaningWater_{WaterMeaning::Disabled};
1024 PressureMeaning primaryVarsMeaningPressure_{PressureMeaning::Po};
1025 GasMeaning primaryVarsMeaningGas_{GasMeaning::Disabled};
1026 BrineMeaning primaryVarsMeaningBrine_{BrineMeaning::Disabled};
1027 SolventMeaning primaryVarsMeaningSolvent_{SolventMeaning::Disabled};
1028 unsigned short pvtRegionIdx_;
1029 Scalar pcFactor_;
1030 static inline Scalar pressureScale_ = 1.0;
1031};
1032
1033} // namespace Opm
1034
1035#endif
Contains the classes required to extend the black-oil model by brine.
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,...
Declares the properties required by the black oil model.
Contains the classes required to extend the black-oil model by solvents.
Contains the high level supplements required to extend the black oil model by brine.
Definition: blackoilbrinemodules.hh:56
static Scalar saltSol(unsigned regionIdx)
Definition: blackoilbrinemodules.hh:333
static const TabulatedFunction & pcfactTable(unsigned satnumRegionIdx)
Definition: blackoilbrinemodules.hh:277
static bool hasPcfactTables()
Definition: blackoilbrinemodules.hh:323
Contains the high level supplements required to extend the black oil model by energy.
Definition: blackoilenergymodules.hh:58
static void assignPrimaryVars(PrimaryVariables &priVars, const FluidState &fluidState)
Assign the energy specific primary variables to a PrimaryVariables object.
Definition: blackoilenergymodules.hh:263
Contains the high level supplements required to extend the black oil model.
Definition: blackoilextbomodules.hh:62
static Value rs(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:354
static Value rv(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:358
Represents the primary variables used by the black-oil model.
Definition: blackoilprimaryvariables.hh:67
bool operator==(const BlackOilPrimaryVariables &rhs) const
Definition: blackoilprimaryvariables.hh:928
WaterMeaning
Definition: blackoilprimaryvariables.hh:127
unsigned pvtRegionIndex() const
Return the index of the region which should be used for PVT properties.
Definition: blackoilprimaryvariables.hh:233
BlackOilPrimaryVariables & operator=(const BlackOilPrimaryVariables &other)=default
SolventMeaning
Definition: blackoilprimaryvariables.hh:153
static void registerParameters()
Definition: blackoilprimaryvariables.hh:192
GasMeaning primaryVarsMeaningGas() const
Return the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:268
PressureMeaning primaryVarsMeaningPressure() const
Return the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:254
void serializeOp(Serializer &serializer)
Definition: blackoilprimaryvariables.hh:916
static void init()
Definition: blackoilprimaryvariables.hh:186
void setPrimaryVarsMeaningSolvent(SolventMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:295
BlackOilPrimaryVariables()
Definition: blackoilprimaryvariables.hh:159
static BlackOilPrimaryVariables serializationTestObject()
Definition: blackoilprimaryvariables.hh:170
SolventMeaning primaryVarsMeaningSolvent() const
Definition: blackoilprimaryvariables.hh:288
bool chopAndNormalizeSaturations()
Definition: blackoilprimaryvariables.hh:841
bool adaptPrimaryVariables(const Problem &problem, unsigned globalDofIdx, Scalar swMaximum, Scalar thresholdWaterFilledCell, Scalar eps=0.0)
Adapt the interpretation of the switching variables to be physically meaningful.
Definition: blackoilprimaryvariables.hh:473
void setPrimaryVarsMeaningPressure(PressureMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:261
WaterMeaning primaryVarsMeaningWater() const
Return the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:240
void setPrimaryVarsMeaningWater(WaterMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:247
void assignNaive(const FluidState &fluidState)
Directly retrieve the primary variables from an arbitrary fluid state.
Definition: blackoilprimaryvariables.hh:302
Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx, LinearizationType linearizationType=LinearizationType()) const
Definition: blackoilprimaryvariables.hh:202
BrineMeaning
Definition: blackoilprimaryvariables.hh:147
void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:285
BlackOilPrimaryVariables(const BlackOilPrimaryVariables &value)=default
Copy constructor.
BrineMeaning primaryVarsMeaningBrine() const
Definition: blackoilprimaryvariables.hh:278
void setPressureScale(Scalar val)
Definition: blackoilprimaryvariables.hh:198
void setPrimaryVarsMeaningGas(GasMeaning newMeaning)
Set the interpretation which should be applied to the switching primary variables.
Definition: blackoilprimaryvariables.hh:275
PressureMeaning
Definition: blackoilprimaryvariables.hh:134
void setPvtRegionIndex(unsigned value)
Set the index of the region which should be used for PVT properties.
Definition: blackoilprimaryvariables.hh:227
void checkDefined() const
< Import base class assignment operators.
Definition: blackoilprimaryvariables.hh:896
GasMeaning
Definition: blackoilprimaryvariables.hh:140
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilsolventmodules.hh:68
static bool isSolubleInWater()
Definition: blackoilsolventmodules.hh:523
static Value solubilityLimit(unsigned pvtIdx, const Value &temperature, const Value &pressure, const Value &saltConcentration)
Definition: blackoilsolventmodules.hh:505
Represents the primary variables used by the a model.
Definition: fvbaseprimaryvariables.hh:53
SimpleModularFluidState< double, 3, 3, FluidSystemSimple, false, false, false, false, true, false, false, false > SatOnlyFluidState
Definition: EquilibrationHelpers_impl.hpp:53
Definition: blackoilnewtonmethodparams.hpp:31
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
Definition: linearizationtype.hh:34
Definition: blackoilprimaryvariables.hh:54
static constexpr Scalar value
Definition: blackoilprimaryvariables.hh:54