24#include <opm/common/TimingMacros.hpp>
26#include <opm/common/utility/numeric/RootFinders.hpp>
28#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
29#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
30#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
34#include <fmt/format.h>
55namespace Miscibility {
57template<
class Flu
idSystem>
59 const std::vector<Scalar>& depth,
60 const std::vector<Scalar>& rs)
61 : pvtRegionIdx_(pvtRegionIdx)
62 , rsVsDepth_(depth, rs)
66template<
class Flu
idSystem>
74 const auto sat_rs = satRs(press, temp);
75 if (satGas > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
79 if (rsVsDepth_.xMin() > depth)
80 return std::min(sat_rs, rsVsDepth_.valueAt(0));
81 else if (rsVsDepth_.xMax() < depth)
82 return std::min(sat_rs, rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1));
83 return std::min(sat_rs, rsVsDepth_.eval(depth,
false));
87template<
class Flu
idSystem>
90satRs(
const Scalar press,
const Scalar temp)
const
92 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
95template<
class Flu
idSystem>
97 const std::vector<Scalar>& depth,
98 const std::vector<Scalar>& pbub)
99 : pvtRegionIdx_(pvtRegionIdx)
100 , pbubVsDepth_(depth, pbub)
104template<
class Flu
idSystem>
110 const Scalar satGas)
const
114 if (pbubVsDepth_.xMin() > depth)
115 press = pbubVsDepth_.valueAt(0);
116 else if (pbubVsDepth_.xMax() < depth)
117 press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1);
119 press = pbubVsDepth_.eval(depth,
false);
121 return satRs(std::min(press, cellPress), temp);
124template<
class Flu
idSystem>
127satRs(
const Scalar press,
const Scalar temp)
const
129 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
132template<
class Flu
idSystem>
134 const std::vector<Scalar>& depth,
135 const std::vector<Scalar>& pdew)
136 : pvtRegionIdx_(pvtRegionIdx)
137 , pdewVsDepth_(depth, pdew)
141template<
class Flu
idSystem>
142typename PDVD<FluidSystem>::Scalar
145 const Scalar cellPress,
147 const Scalar satOil)
const
149 Scalar press = cellPress;
151 if (pdewVsDepth_.xMin() > depth)
152 press = pdewVsDepth_.valueAt(0);
153 else if (pdewVsDepth_.xMax() < depth)
154 press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1);
156 press = pdewVsDepth_.eval(depth,
false);
158 return satRv(std::min(press, cellPress), temp);
161template<
class Flu
idSystem>
162typename PDVD<FluidSystem>::Scalar
164satRv(
const Scalar press,
const Scalar temp)
const
166 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
169template<
class Flu
idSystem>
171 const std::vector<Scalar>& depth,
172 const std::vector<Scalar>& rv)
173 : pvtRegionIdx_(pvtRegionIdx)
174 , rvVsDepth_(depth, rv)
178template<
class Flu
idSystem>
179typename RvVD<FluidSystem>::Scalar
184 const Scalar satOil)
const
186 if (satOil < - std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
187 throw std::logic_error {
188 "Must not pass negative oil saturation"
191 const auto sat_rv = satRv(press, temp);
192 if (satOil > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
196 if (rvVsDepth_.xMin() > depth)
197 return std::min(sat_rv, rvVsDepth_.valueAt(0));
198 else if (rvVsDepth_.xMax() < depth)
199 return std::min(sat_rv, rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1));
200 return std::min(sat_rv, rvVsDepth_.eval(depth,
false));
204template<
class Flu
idSystem>
205typename RvVD<FluidSystem>::Scalar
207satRv(
const Scalar press,
const Scalar temp)
const
209 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
212template<
class Flu
idSystem>
214 const std::vector<Scalar>& depth,
215 const std::vector<Scalar>& rvw)
216 : pvtRegionIdx_(pvtRegionIdx)
217 , rvwVsDepth_(depth, rvw)
221template<
class Flu
idSystem>
222typename RvwVD<FluidSystem>::Scalar
227 const Scalar satWat)
const
229 if (satWat < - std::sqrt(std::numeric_limits<double>::epsilon())) {
230 throw std::logic_error {
231 "Must not pass negative water saturation"
235 const auto sat_rvw = satRvw(press, temp);
236 if (satWat > std::sqrt(std::numeric_limits<double>::epsilon())) {
240 if (rvwVsDepth_.xMin() > depth)
241 return std::min(sat_rvw,rvwVsDepth_.valueAt(0));
242 else if (rvwVsDepth_.xMax() < depth)
243 return std::min(sat_rvw, rvwVsDepth_.valueAt(rvwVsDepth_.numSamples() - 1));
244 return std::min(sat_rvw, rvwVsDepth_.eval(depth,
false));
248template<
class Flu
idSystem>
249typename RvwVD<FluidSystem>::Scalar
251satRvw(
const Scalar press,
const Scalar temp)
const
253 return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);
256template<
class Flu
idSystem>
258RsSatAtContact(
const int pvtRegionIdx,
const Scalar pContact,
const Scalar T_contact)
259 : pvtRegionIdx_(pvtRegionIdx)
261 rsSatContact_ = satRs(pContact, T_contact);
264template<
class Flu
idSystem>
265typename RsSatAtContact<FluidSystem>::Scalar
270 const Scalar satGas)
const
272 if (satGas > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
273 return satRs(press, temp);
276 return std::min(satRs(press, temp), rsSatContact_);
280template<
class Flu
idSystem>
281typename RsSatAtContact<FluidSystem>::Scalar
283satRs(
const Scalar press,
const Scalar temp)
const
285 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
288template<
class Flu
idSystem>
290RvSatAtContact(
const int pvtRegionIdx,
const Scalar pContact,
const Scalar T_contact)
291 : pvtRegionIdx_(pvtRegionIdx)
293 rvSatContact_ = satRv(pContact, T_contact);
296template<
class Flu
idSystem>
297typename RvSatAtContact<FluidSystem>::Scalar
302 const Scalar satOil)
const
304 if (satOil > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
305 return satRv(press, temp);
308 return std::min(satRv(press, temp), rvSatContact_);
312template<
class Flu
idSystem>
313typename RvSatAtContact<FluidSystem>::Scalar
315satRv(
const Scalar press,
const Scalar temp)
const
317 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
320template<
class Flu
idSystem>
322RvwSatAtContact(
const int pvtRegionIdx,
const Scalar pContact,
const Scalar T_contact)
323 : pvtRegionIdx_(pvtRegionIdx)
325 rvwSatContact_ = satRvw(pContact, T_contact);
328template<
class Flu
idSystem>
329typename RvwSatAtContact<FluidSystem>::Scalar
334 const Scalar satWat)
const
336 if (satWat > std::sqrt(std::numeric_limits<Scalar>::epsilon())) {
337 return satRvw(press, temp);
340 return std::min(satRvw(press, temp), rvwSatContact_);
344template<
class Flu
idSystem>
345typename RvwSatAtContact<FluidSystem>::Scalar
347satRvw(
const Scalar press,
const Scalar temp)
const
349 return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);;
354template<
class Scalar>
359 const TabulatedFunction& tempVdTable,
360 const TabulatedFunction& saltVdTable,
366 , tempVdTable_ (tempVdTable)
367 , saltVdTable_ (saltVdTable)
372template<
class Scalar>
375 return this->rec_.datumDepth();
378template<
class Scalar>
381 return this->rec_.datumDepthPressure();
384template<
class Scalar>
387 return this->rec_.waterOilContactDepth();
390template<
class Scalar>
393 return this->rec_.waterOilContactCapillaryPressure();
396template<
class Scalar>
399 return this->rec_.gasOilContactDepth();
402template<
class Scalar>
405 return this->rec_.gasOilContactCapillaryPressure();
408template<
class Scalar>
411 return this->rec_.initializationTargetAccuracy();
414template<
class Scalar>
421template<
class Scalar>
428template<
class Scalar>
435template<
class Scalar>
436const typename EquilReg<Scalar>::TabulatedFunction&
442template<
class Scalar>
443const typename EquilReg<Scalar>::TabulatedFunction&
449template<
class Scalar>
452 return this->pvtIdx_;
455template<
class Flu
idSystem,
class MaterialLawManager>
457PcEq(
const MaterialLawManager& materialLawManager,
461 : materialLawManager_(materialLawManager)
464 , targetPc_(targetPc)
468template<
class Flu
idSystem,
class MaterialLawManager>
473 const auto& matParams = materialLawManager_.materialLawParams(cell_);
475 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
476 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
477 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
478 fluidState.setSaturation(phase_, s);
481 assert(phase_ != FluidSystem::oilPhaseIdx);
483 if (phase_ == FluidSystem::gasPhaseIdx) {
484 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0 - s);
485 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - s);
487 std::array<Scalar, FluidSystem::numPhases> pc{0.0};
488 using MaterialLaw =
typename MaterialLawManager::MaterialLaw;
489 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
490 Scalar sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
491 Scalar pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
492 return pcPhase - targetPc_;
495template<
class Flu
idSystem,
class MaterialLawManager>
497PcEqSum(
const MaterialLawManager& materialLawManager,
502 : materialLawManager_(materialLawManager)
506 , targetPc_(targetPc)
510template<
class Flu
idSystem,
class MaterialLawManager>
515 const auto& matParams = materialLawManager_.materialLawParams(cell_);
517 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
518 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
519 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
520 fluidState.setSaturation(phase1_, s);
521 fluidState.setSaturation(phase2_, 1.0 - s);
523 std::array<Scalar, FluidSystem::numPhases> pc {0.0};
525 using MaterialLaw =
typename MaterialLawManager::MaterialLaw;
526 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
527 Scalar sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
528 Scalar pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
529 Scalar sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
530 Scalar pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
531 return pc1 + pc2 - targetPc_;
534template <
class Flu
idSystem,
class MaterialLawManager>
535typename FluidSystem::Scalar
537 const int phase,
const int cell)
539 const auto& scaledDrainageInfo =
540 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
544 case FluidSystem::waterPhaseIdx:
545 return scaledDrainageInfo.Swl;
547 case FluidSystem::gasPhaseIdx:
548 return scaledDrainageInfo.Sgl;
550 case FluidSystem::oilPhaseIdx:
551 throw std::runtime_error(
"Min saturation not implemented for oil phase.");
554 throw std::runtime_error(
"Unknown phaseIdx .");
559template <
class Flu
idSystem,
class MaterialLawManager>
560typename FluidSystem::Scalar
562 const int phase,
const int cell)
564 const auto& scaledDrainageInfo =
565 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
569 case FluidSystem::waterPhaseIdx:
570 return scaledDrainageInfo.Swu;
572 case FluidSystem::gasPhaseIdx:
573 return scaledDrainageInfo.Sgu;
575 case FluidSystem::oilPhaseIdx:
576 throw std::runtime_error(
"Max saturation not implemented for oil phase.");
579 throw std::runtime_error(
"Unknown phaseIdx .");
584template <
class Flu
idSystem,
class MaterialLawManager>
585typename FluidSystem::Scalar
589 const typename FluidSystem::Scalar targetPc,
590 const bool increasing)
592 using Scalar =
typename FluidSystem::Scalar;
594 Scalar s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
595 Scalar s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
601 if (!std::isfinite(f0 + f1))
602 throw std::logic_error(fmt::format(
"The capillary pressure values {} and {} are not finite", f0, f1));
609 const Scalar tol = 1e-10;
611 const int maxIter = -2*
static_cast<int>(std::log2(tol)) + 10;
612 int usedIterations = -1;
613 const Scalar root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
617template<
class Flu
idSystem,
class MaterialLawManager>
618typename FluidSystem::Scalar
623 const typename FluidSystem::Scalar targetPc)
625 using Scalar =
typename FluidSystem::Scalar;
627 Scalar s0 = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
628 Scalar s1 = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
639 assert(f0 > 0.0 && f1 < 0.0);
640 const Scalar tol = 1e-10;
642 const int maxIter = -2*
static_cast<int>(std::log2(tol)) + 10;
643 int usedIterations = -1;
644 const Scalar root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
648template<
class Flu
idSystem,
class MaterialLawManager>
649typename FluidSystem::Scalar
651 const typename FluidSystem::Scalar cellDepth,
652 const typename FluidSystem::Scalar contactDepth,
655 const bool increasing)
657 using Scalar =
typename FluidSystem::Scalar;
658 const Scalar s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
659 const Scalar s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
661 if (cellDepth < contactDepth) {
669template<
class Flu
idSystem,
class MaterialLawManager>
670bool isConstPc(
const MaterialLawManager& materialLawManager,
674 using Scalar =
typename FluidSystem::Scalar;
677 const Scalar f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
678 const Scalar f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
679 return std::abs(f0 - f1) < std::numeric_limits<Scalar>::epsilon();
Auxiliary routines that to solve the ODEs that emerge from the hydrostatic equilibrium problem.
Scalar datum() const
Definition: EquilibrationHelpers_impl.hpp:373
int pvtIdx() const
Definition: EquilibrationHelpers_impl.hpp:450
EquilReg(const EquilRecord &rec, std::shared_ptr< Miscibility::RsFunction< Scalar > > rs, std::shared_ptr< Miscibility::RsFunction< Scalar > > rv, std::shared_ptr< Miscibility::RsFunction< Scalar > > rvw, const TabulatedFunction &tempVdTable, const TabulatedFunction &saltVdTable, const int pvtIdx)
Definition: EquilibrationHelpers_impl.hpp:355
Scalar pcgoGoc() const
Definition: EquilibrationHelpers_impl.hpp:403
int equilibrationAccuracy() const
Definition: EquilibrationHelpers_impl.hpp:409
const CalcEvaporation & evaporationCalculator() const
Definition: EquilibrationHelpers_impl.hpp:423
Scalar zgoc() const
Definition: EquilibrationHelpers_impl.hpp:397
Scalar pressure() const
Definition: EquilibrationHelpers_impl.hpp:379
const TabulatedFunction & saltVdTable() const
Definition: EquilibrationHelpers_impl.hpp:437
Scalar pcowWoc() const
Definition: EquilibrationHelpers_impl.hpp:391
const TabulatedFunction & tempVdTable() const
Definition: EquilibrationHelpers_impl.hpp:444
Scalar zwoc() const
Definition: EquilibrationHelpers_impl.hpp:385
const CalcDissolution & dissolutionCalculator() const
Definition: EquilibrationHelpers_impl.hpp:416
const CalcWaterEvaporation & waterEvaporationCalculator() const
Definition: EquilibrationHelpers_impl.hpp:430
Definition: EquilibrationHelpers.hpp:223
Scalar operator()(const Scalar depth, const Scalar cellPress, const Scalar temp, const Scalar satGas=0.0) const override
Definition: EquilibrationHelpers_impl.hpp:107
PBVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &pbub)
Definition: EquilibrationHelpers_impl.hpp:96
typename FluidSystem::Scalar Scalar
Definition: EquilibrationHelpers.hpp:225
Definition: EquilibrationHelpers.hpp:275
Scalar operator()(const Scalar depth, const Scalar cellPress, const Scalar temp, const Scalar satOil=0.0) const override
Definition: EquilibrationHelpers_impl.hpp:144
PDVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &pdew)
Definition: EquilibrationHelpers_impl.hpp:133
Definition: EquilibrationHelpers.hpp:101
Definition: EquilibrationHelpers.hpp:170
Scalar operator()(const Scalar depth, const Scalar press, const Scalar temp, const Scalar satGas=0.0) const override
Definition: EquilibrationHelpers_impl.hpp:69
typename FluidSystem::Scalar Scalar
Definition: EquilibrationHelpers.hpp:172
RsVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &rs)
Definition: EquilibrationHelpers_impl.hpp:58
Definition: EquilibrationHelpers.hpp:327
Scalar operator()(const Scalar depth, const Scalar press, const Scalar temp, const Scalar satOil=0.0) const override
Definition: EquilibrationHelpers_impl.hpp:181
RvVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &rv)
Definition: EquilibrationHelpers_impl.hpp:170
Definition: EquilibrationHelpers.hpp:378
Scalar operator()(const Scalar depth, const Scalar press, const Scalar temp, const Scalar satWat=0.0) const override
Definition: EquilibrationHelpers_impl.hpp:224
RvwVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &rvw)
Definition: EquilibrationHelpers_impl.hpp:213
BlackOilFluidSystem< double > FluidSystemSimple
Definition: EquilibrationHelpers_impl.hpp:39
FluidSystem::Scalar satFromDepth(const MaterialLawManager &materialLawManager, const typename FluidSystem::Scalar cellDepth, const typename FluidSystem::Scalar contactDepth, const int phase, const int cell, const bool increasing=false)
Compute saturation from depth. Used for constant capillary pressure function.
Definition: EquilibrationHelpers_impl.hpp:650
FluidSystem::Scalar satFromPc(const MaterialLawManager &materialLawManager, const int phase, const int cell, const typename FluidSystem::Scalar targetPc, const bool increasing=false)
Definition: EquilibrationHelpers_impl.hpp:586
SimpleModularFluidState< double, 3, 3, FluidSystemSimple, false, false, false, false, true, false, false, false > SatOnlyFluidState
Definition: EquilibrationHelpers_impl.hpp:53
FluidSystem::Scalar minSaturations(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Definition: EquilibrationHelpers_impl.hpp:536
FluidSystem::Scalar maxSaturations(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Definition: EquilibrationHelpers_impl.hpp:561
bool isConstPc(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition: EquilibrationHelpers_impl.hpp:670
FluidSystem::Scalar satFromSumOfPcs(const MaterialLawManager &materialLawManager, const int phase1, const int phase2, const int cell, const typename FluidSystem::Scalar targetPc)
Definition: EquilibrationHelpers_impl.hpp:619
Definition: blackoilboundaryratevector.hh:37
Definition: EquilibrationHelpers.hpp:740
Scalar operator()(Scalar s) const
Definition: EquilibrationHelpers_impl.hpp:471
typename FluidSystem::Scalar Scalar
Definition: EquilibrationHelpers.hpp:741
PcEq(const MaterialLawManager &materialLawManager, const int phase, const int cell, const Scalar targetPc)
Definition: EquilibrationHelpers_impl.hpp:457
Definition: EquilibrationHelpers.hpp:781
PcEqSum(const MaterialLawManager &materialLawManager, const int phase1, const int phase2, const int cell, const Scalar targetPc)
Definition: EquilibrationHelpers_impl.hpp:497
Scalar operator()(Scalar s) const
Definition: EquilibrationHelpers_impl.hpp:513
typename FluidSystem::Scalar Scalar
Definition: EquilibrationHelpers.hpp:782