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<double>& depth,
60 const std::vector<double>& rs)
61 : pvtRegionIdx_(pvtRegionIdx)
62 , rsVsDepth_(depth, rs)
66template<
class Flu
idSystem>
71 const double satGas)
const
73 const auto sat_rs = satRs(press, temp);
74 if (satGas > std::sqrt(std::numeric_limits<double>::epsilon())) {
78 if (rsVsDepth_.xMin() > depth)
79 return std::min(sat_rs, rsVsDepth_.valueAt(0));
80 else if (rsVsDepth_.xMax() < depth)
81 return std::min(sat_rs, rsVsDepth_.valueAt(rsVsDepth_.numSamples() - 1));
82 return std::min(sat_rs, rsVsDepth_.eval(depth,
false));
86template<
class Flu
idSystem>
89 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
92template<
class Flu
idSystem>
94 const std::vector<double>& depth,
95 const std::vector<double>& pbub)
96 : pvtRegionIdx_(pvtRegionIdx)
97 , pbubVsDepth_(depth, pbub)
101template<
class Flu
idSystem>
104 const double cellPress,
106 const double satGas)
const
108 double press = cellPress;
110 if (pbubVsDepth_.xMin() > depth)
111 press = pbubVsDepth_.valueAt(0);
112 else if (pbubVsDepth_.xMax() < depth)
113 press = pbubVsDepth_.valueAt(pbubVsDepth_.numSamples() - 1);
115 press = pbubVsDepth_.eval(depth,
false);
117 return satRs(std::min(press, cellPress), temp);
120template<
class Flu
idSystem>
122satRs(
const double press,
const double temp)
const
124 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
127template<
class Flu
idSystem>
129 const std::vector<double>& depth,
130 const std::vector<double>& pdew)
131 : pvtRegionIdx_(pvtRegionIdx)
132 , pdewVsDepth_(depth, pdew)
136template<
class Flu
idSystem>
139 const double cellPress,
141 const double satOil)
const
143 double press = cellPress;
145 if (pdewVsDepth_.xMin() > depth)
146 press = pdewVsDepth_.valueAt(0);
147 else if (pdewVsDepth_.xMax() < depth)
148 press = pdewVsDepth_.valueAt(pdewVsDepth_.numSamples() - 1);
150 press = pdewVsDepth_.eval(depth,
false);
152 return satRv(std::min(press, cellPress), temp);
155template<
class Flu
idSystem>
157satRv(
const double press,
const double temp)
const
159 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
163template<
class Flu
idSystem>
165 const std::vector<double>& depth,
166 const std::vector<double>& rv)
167 : pvtRegionIdx_(pvtRegionIdx)
168 , rvVsDepth_(depth, rv)
172template<
class Flu
idSystem>
177 const double satOil)
const
179 if (satOil < - std::sqrt(std::numeric_limits<double>::epsilon())) {
180 throw std::logic_error {
181 "Must not pass negative oil saturation"
184 const auto sat_rv = satRv(press, temp);
185 if (satOil > std::sqrt(std::numeric_limits<double>::epsilon())) {
189 if (rvVsDepth_.xMin() > depth)
190 return std::min(sat_rv, rvVsDepth_.valueAt(0));
191 else if (rvVsDepth_.xMax() < depth)
192 return std::min(sat_rv, rvVsDepth_.valueAt(rvVsDepth_.numSamples() - 1));
193 return std::min(sat_rv, rvVsDepth_.eval(depth,
false));
197template<
class Flu
idSystem>
199satRv(
const double press,
const double temp)
const
201 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
205template<
class Flu
idSystem>
207 const std::vector<double>& depth,
208 const std::vector<double>& rvw)
209 : pvtRegionIdx_(pvtRegionIdx)
210 , rvwVsDepth_(depth, rvw)
214template<
class Flu
idSystem>
219 const double satWat)
const
221 if (satWat < - std::sqrt(std::numeric_limits<double>::epsilon())) {
222 throw std::logic_error {
223 "Must not pass negative water saturation"
227 const auto sat_rvw = satRvw(press, temp);
228 if (satWat > std::sqrt(std::numeric_limits<double>::epsilon())) {
232 if (rvwVsDepth_.xMin() > depth)
233 return std::min(sat_rvw,rvwVsDepth_.valueAt(0));
234 else if (rvwVsDepth_.xMax() < depth)
235 return std::min(sat_rvw, rvwVsDepth_.valueAt(rvwVsDepth_.numSamples() - 1));
236 return std::min(sat_rvw, rvwVsDepth_.eval(depth,
false));
240template<
class Flu
idSystem>
242satRvw(
const double press,
const double temp)
const
244 return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);
248template<
class Flu
idSystem>
250RsSatAtContact(
const int pvtRegionIdx,
const double pContact,
const double T_contact)
251 : pvtRegionIdx_(pvtRegionIdx)
253 rsSatContact_ = satRs(pContact, T_contact);
256template<
class Flu
idSystem>
261 const double satGas)
const
263 if (satGas > std::sqrt(std::numeric_limits<double>::epsilon())) {
264 return satRs(press, temp);
267 return std::min(satRs(press, temp), rsSatContact_);
271template<
class Flu
idSystem>
273satRs(
const double press,
const double temp)
const
275 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
278template<
class Flu
idSystem>
280RvSatAtContact(
const int pvtRegionIdx,
const double pContact,
const double T_contact)
281 : pvtRegionIdx_(pvtRegionIdx)
283 rvSatContact_ = satRv(pContact, T_contact);
286template<
class Flu
idSystem>
291 const double satOil)
const
293 if (satOil > std::sqrt(std::numeric_limits<double>::epsilon())) {
294 return satRv(press, temp);
297 return std::min(satRv(press, temp), rvSatContact_);
301template<
class Flu
idSystem>
303satRv(
const double press,
const double temp)
const
305 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
308template<
class Flu
idSystem>
310RvwSatAtContact(
const int pvtRegionIdx,
const double pContact,
const double T_contact)
311 : pvtRegionIdx_(pvtRegionIdx)
313 rvwSatContact_ = satRvw(pContact, T_contact);
316template<
class Flu
idSystem>
321 const double satWat)
const
323 if (satWat > std::sqrt(std::numeric_limits<double>::epsilon())) {
324 return satRvw(press, temp);
327 return std::min(satRvw(press, temp), rvwSatContact_);
331template<
class Flu
idSystem>
333satRvw(
const double press,
const double temp)
const
335 return FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_, temp, press);;
341 std::shared_ptr<Miscibility::RsFunction> rs,
342 std::shared_ptr<Miscibility::RsFunction> rv,
343 std::shared_ptr<Miscibility::RsFunction> rvw,
344 const TabulatedFunction& tempVdTable,
345 const TabulatedFunction& saltVdTable,
351 , tempVdTable_ (tempVdTable)
352 , saltVdTable_ (saltVdTable)
359 return this->rec_.datumDepth();
364 return this->rec_.datumDepthPressure();
369 return this->rec_.waterOilContactDepth();
374 return this->rec_.waterOilContactCapillaryPressure();
379 return this->rec_.gasOilContactDepth();
384 return this->rec_.gasOilContactCapillaryPressure();
389 return this->rec_.initializationTargetAccuracy();
410const EquilReg::TabulatedFunction&
416const EquilReg::TabulatedFunction&
424 return this->pvtIdx_;
427template<
class Flu
idSystem,
class MaterialLawManager>
429PcEq(
const MaterialLawManager& materialLawManager,
432 const double targetPc)
433 : materialLawManager_(materialLawManager),
440template<
class Flu
idSystem,
class MaterialLawManager>
444 const auto& matParams = materialLawManager_.materialLawParams(cell_);
446 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
447 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
448 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
449 fluidState.setSaturation(phase_, s);
452 assert(phase_ != FluidSystem::oilPhaseIdx);
454 if (phase_ == FluidSystem::gasPhaseIdx) {
455 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0 - s);
456 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - s);
458 std::array<double, FluidSystem::numPhases> pc{0.0};
459 using MaterialLaw =
typename MaterialLawManager::MaterialLaw;
460 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
461 double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
462 double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
463 return pcPhase - targetPc_;
466template<
class Flu
idSystem,
class MaterialLawManager>
468PcEqSum(
const MaterialLawManager& materialLawManager,
472 const double targetPc)
473 : materialLawManager_(materialLawManager),
481template<
class Flu
idSystem,
class MaterialLawManager>
485 const auto& matParams = materialLawManager_.materialLawParams(cell_);
487 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
488 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
489 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
490 fluidState.setSaturation(phase1_, s);
491 fluidState.setSaturation(phase2_, 1.0 - s);
493 std::array<double, FluidSystem::numPhases> pc {0.0};
495 using MaterialLaw =
typename MaterialLawManager::MaterialLaw;
496 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
497 double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
498 double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
499 double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
500 double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
501 return pc1 + pc2 - targetPc_;
504template <
class Flu
idSystem,
class MaterialLawManager>
506 const int phase,
const int cell)
508 const auto& scaledDrainageInfo =
509 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
513 case FluidSystem::waterPhaseIdx:
514 return scaledDrainageInfo.Swl;
516 case FluidSystem::gasPhaseIdx:
517 return scaledDrainageInfo.Sgl;
519 case FluidSystem::oilPhaseIdx:
520 throw std::runtime_error(
"Min saturation not implemented for oil phase.");
523 throw std::runtime_error(
"Unknown phaseIdx .");
528template <
class Flu
idSystem,
class MaterialLawManager>
530 const int phase,
const int cell)
532 const auto& scaledDrainageInfo =
533 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
537 case FluidSystem::waterPhaseIdx:
538 return scaledDrainageInfo.Swu;
540 case FluidSystem::gasPhaseIdx:
541 return scaledDrainageInfo.Sgu;
543 case FluidSystem::oilPhaseIdx:
544 throw std::runtime_error(
"Max saturation not implemented for oil phase.");
547 throw std::runtime_error(
"Unknown phaseIdx .");
552template <
class Flu
idSystem,
class MaterialLawManager>
553double satFromPc(
const MaterialLawManager& materialLawManager,
556 const double targetPc,
557 const bool increasing)
560 double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
561 double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
567 if (!std::isfinite(f0 + f1))
568 throw std::logic_error(fmt::format(
"The capillary pressure values {} and {} are not finite", f0, f1));
575 const double tol = 1e-10;
577 const int maxIter = -2*
static_cast<int>(std::log2(tol)) + 10;
578 int usedIterations = -1;
579 const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
583template<
class Flu
idSystem,
class MaterialLawManager>
588 const double targetPc)
591 double s0 = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
592 double s1 = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
603 assert(f0 > 0.0 && f1 < 0.0);
604 const double tol = 1e-10;
606 const int maxIter = -2*
static_cast<int>(std::log2(tol)) + 10;
607 int usedIterations = -1;
608 const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
612template<
class Flu
idSystem,
class MaterialLawManager>
614 const double cellDepth,
615 const double contactDepth,
618 const bool increasing)
620 const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
621 const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
623 if (cellDepth < contactDepth) {
631template<
class Flu
idSystem,
class MaterialLawManager>
632bool isConstPc(
const MaterialLawManager& materialLawManager,
638 const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
639 const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
640 return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
Auxiliary routines that to solve the ODEs that emerge from the hydrostatic equilibrium problem.
const CalcDissolution & dissolutionCalculator() const
Definition: EquilibrationHelpers_impl.hpp:393
double datum() const
Definition: EquilibrationHelpers_impl.hpp:357
double pcgoGoc() const
Definition: EquilibrationHelpers_impl.hpp:382
int equilibrationAccuracy() const
Definition: EquilibrationHelpers_impl.hpp:387
double zgoc() const
Definition: EquilibrationHelpers_impl.hpp:377
const CalcWaterEvaporation & waterEvaporationCalculator() const
Definition: EquilibrationHelpers_impl.hpp:405
const TabulatedFunction & tempVdTable() const
Definition: EquilibrationHelpers_impl.hpp:417
const CalcEvaporation & evaporationCalculator() const
Definition: EquilibrationHelpers_impl.hpp:399
double pcowWoc() const
Definition: EquilibrationHelpers_impl.hpp:372
double zwoc() const
Definition: EquilibrationHelpers_impl.hpp:367
double pressure() const
Definition: EquilibrationHelpers_impl.hpp:362
int pvtIdx() const
Definition: EquilibrationHelpers_impl.hpp:422
EquilReg(const EquilRecord &rec, std::shared_ptr< Miscibility::RsFunction > rs, std::shared_ptr< Miscibility::RsFunction > rv, std::shared_ptr< Miscibility::RsFunction > rvw, const TabulatedFunction &tempVdTable, const TabulatedFunction &saltVdTable, const int pvtIdx)
Definition: EquilibrationHelpers_impl.hpp:340
const TabulatedFunction & saltVdTable() const
Definition: EquilibrationHelpers_impl.hpp:411
Definition: EquilibrationHelpers.hpp:220
double operator()(const double depth, const double cellPress, const double temp, const double satGas=0.0) const
Definition: EquilibrationHelpers_impl.hpp:103
PBVD(const int pvtRegionIdx, const std::vector< double > &depth, const std::vector< double > &pbub)
Definition: EquilibrationHelpers_impl.hpp:93
Definition: EquilibrationHelpers.hpp:271
PDVD(const int pvtRegionIdx, const std::vector< double > &depth, const std::vector< double > &pdew)
Definition: EquilibrationHelpers_impl.hpp:128
double operator()(const double depth, const double cellPress, const double temp, const double satOil=0.0) const
Definition: EquilibrationHelpers_impl.hpp:138
Definition: EquilibrationHelpers.hpp:100
Definition: EquilibrationHelpers.hpp:168
RsVD(const int pvtRegionIdx, const std::vector< double > &depth, const std::vector< double > &rs)
Definition: EquilibrationHelpers_impl.hpp:58
double operator()(const double depth, const double press, const double temp, const double satGas=0.0) const
Definition: EquilibrationHelpers_impl.hpp:68
Definition: EquilibrationHelpers.hpp:322
double operator()(const double depth, const double press, const double temp, const double satOil=0.0) const
Definition: EquilibrationHelpers_impl.hpp:174
RvVD(const int pvtRegionIdx, const std::vector< double > &depth, const std::vector< double > &rv)
Definition: EquilibrationHelpers_impl.hpp:164
Definition: EquilibrationHelpers.hpp:372
double operator()(const double depth, const double press, const double temp, const double satWat=0.0) const
Definition: EquilibrationHelpers_impl.hpp:216
RvwVD(const int pvtRegionIdx, const std::vector< double > &depth, const std::vector< double > &rvw)
Definition: EquilibrationHelpers_impl.hpp:206
BlackOilFluidSystem< double > FluidSystemSimple
Definition: EquilibrationHelpers_impl.hpp:39
double satFromDepth(const MaterialLawManager &materialLawManager, const double cellDepth, const double 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:613
double satFromSumOfPcs(const MaterialLawManager &materialLawManager, const int phase1, const int phase2, const int cell, const double targetPc)
Definition: EquilibrationHelpers_impl.hpp:584
double satFromPc(const MaterialLawManager &materialLawManager, const int phase, const int cell, const double targetPc, const bool increasing=false)
Definition: EquilibrationHelpers_impl.hpp:553
SimpleModularFluidState< double, 3, 3, FluidSystemSimple, false, false, false, false, true, false, false, false > SatOnlyFluidState
Definition: EquilibrationHelpers_impl.hpp:53
bool isConstPc(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition: EquilibrationHelpers_impl.hpp:632
double maxSaturations(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Definition: EquilibrationHelpers_impl.hpp:529
double minSaturations(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Definition: EquilibrationHelpers_impl.hpp:505
Definition: BlackoilPhases.hpp:27
Definition: EquilibrationHelpers.hpp:723
double operator()(double s) const
Definition: EquilibrationHelpers_impl.hpp:442
PcEq(const MaterialLawManager &materialLawManager, const int phase, const int cell, const double targetPc)
Definition: EquilibrationHelpers_impl.hpp:429
Definition: EquilibrationHelpers.hpp:760
PcEqSum(const MaterialLawManager &materialLawManager, const int phase1, const int phase2, const int cell, const double targetPc)
Definition: EquilibrationHelpers_impl.hpp:468
double operator()(double s) const
Definition: EquilibrationHelpers_impl.hpp:483