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>
67typename RsVD<FluidSystem>::Scalar
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>
88typename RsVD<FluidSystem>::Scalar
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>
105typename PBVD<FluidSystem>::Scalar
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);
124template<
class Flu
idSystem>
125typename PBVD<FluidSystem>::Scalar
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
147 const Scalar
satOil)
const
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);
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"
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"
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>
259 : pvtRegionIdx_(pvtRegionIdx)
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>
291 : pvtRegionIdx_(pvtRegionIdx)
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>
323 : pvtRegionIdx_(pvtRegionIdx)
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&
444EquilReg<Scalar>::tempVdTable()
const
449template<
class Scalar>
452 return this->pvtIdx_;
455template<
class Flu
idSystem,
class MaterialLawManager>
457PcEq(
const MaterialLawManager& materialLawManager,
461 : materialLawManager_(materialLawManager)
468template<
class Flu
idSystem,
class MaterialLawManager>
469typename PcEq<FluidSystem,MaterialLawManager>::Scalar
470PcEq<FluidSystem,MaterialLawManager>::
471operator()(Scalar s)
const
473 const auto&
matParams = materialLawManager_.materialLawParams(cell_);
474 SatOnlyFluidState fluidState;
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;
495template<
class Flu
idSystem,
class MaterialLawManager>
496PcEqSum<FluidSystem,MaterialLawManager>::
497PcEqSum(
const MaterialLawManager& materialLawManager,
502 : materialLawManager_(materialLawManager)
510template<
class Flu
idSystem,
class MaterialLawManager>
511typename PcEqSum<FluidSystem,MaterialLawManager>::Scalar
512PcEqSum<FluidSystem,MaterialLawManager>::
513operator()(Scalar s)
const
515 const auto&
matParams = materialLawManager_.materialLawParams(cell_);
516 SatOnlyFluidState fluidState;
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
536minSaturations(
const MaterialLawManager& materialLawManager,
537 const int phase,
const int cell)
540 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
544 case FluidSystem::waterPhaseIdx:
547 case FluidSystem::gasPhaseIdx:
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
561maxSaturations(
const MaterialLawManager& materialLawManager,
562 const int phase,
const int cell)
565 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
569 case FluidSystem::waterPhaseIdx:
572 case FluidSystem::gasPhaseIdx:
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,
592 using Scalar =
typename FluidSystem::Scalar;
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 = 1
e-10;
611 const int maxIter = -2*
static_cast<int>(std::log2(
tol)) + 10;
617template<
class Flu
idSystem,
class MaterialLawManager>
618typename FluidSystem::Scalar
623 const typename FluidSystem::Scalar
targetPc)
625 using Scalar =
typename FluidSystem::Scalar;
640 const Scalar
tol = 1
e-10;
642 const int maxIter = -2*
static_cast<int>(std::log2(
tol)) + 10;
648template<
class Flu
idSystem,
class MaterialLawManager>
649typename FluidSystem::Scalar
651 const typename FluidSystem::Scalar
cellDepth,
657 using Scalar =
typename FluidSystem::Scalar;
669template<
class Flu
idSystem,
class MaterialLawManager>
670bool isConstPc(
const MaterialLawManager& materialLawManager,
674 using Scalar =
typename FluidSystem::Scalar;
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.
Aggregate information base of an equilibration region.
Definition EquilibrationHelpers.hpp:609
Scalar datum() const
Datum depth in current region.
Definition EquilibrationHelpers_impl.hpp:373
int pvtIdx() const
Retrieve pvtIdx of the region.
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)
Constructor.
Definition EquilibrationHelpers_impl.hpp:355
Scalar pcgoGoc() const
Gas-oil capillary pressure at gas-oil contact.
Definition EquilibrationHelpers_impl.hpp:403
int equilibrationAccuracy() const
Accuracy/strategy for initial fluid-in-place calculation.
Definition EquilibrationHelpers_impl.hpp:409
const CalcEvaporation & evaporationCalculator() const
Retrieve vapourised oil-gas ratio calculator of current region.
Definition EquilibrationHelpers_impl.hpp:423
Scalar zgoc() const
Depth of gas-oil contact.
Definition EquilibrationHelpers_impl.hpp:397
Scalar pressure() const
Pressure at datum depth in current region.
Definition EquilibrationHelpers_impl.hpp:379
Scalar pcowWoc() const
water-oil capillary pressure at water-oil contact.
Definition EquilibrationHelpers_impl.hpp:391
Scalar zwoc() const
Depth of water-oil contact.
Definition EquilibrationHelpers_impl.hpp:385
const CalcDissolution & dissolutionCalculator() const
Retrieve dissolved gas-oil ratio calculator of current region.
Definition EquilibrationHelpers_impl.hpp:416
const CalcWaterEvaporation & waterEvaporationCalculator() const
Retrieve vapourised water-gas ratio calculator of current region.
Definition EquilibrationHelpers_impl.hpp:430
Type that implements "dissolved gas-oil ratio" tabulated as a function of depth policy.
Definition EquilibrationHelpers.hpp:219
Scalar operator()(const Scalar depth, const Scalar cellPress, const Scalar temp, const Scalar satGas=0.0) const override
Function call.
Definition EquilibrationHelpers_impl.hpp:107
PBVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &pbub)
Constructor.
Definition EquilibrationHelpers_impl.hpp:96
Type that implements "vaporized oil-gas ratio" tabulated as a function of depth policy.
Definition EquilibrationHelpers.hpp:269
Scalar operator()(const Scalar depth, const Scalar cellPress, const Scalar temp, const Scalar satOil=0.0) const override
Function call.
Definition EquilibrationHelpers_impl.hpp:144
PDVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &pdew)
Constructor.
Definition EquilibrationHelpers_impl.hpp:133
Base class for phase mixing functions.
Definition EquilibrationHelpers.hpp:101
Type that implements "dissolved gas-oil ratio" tabulated as a function of depth policy.
Definition EquilibrationHelpers.hpp:168
Scalar operator()(const Scalar depth, const Scalar press, const Scalar temp, const Scalar satGas=0.0) const override
Function call.
Definition EquilibrationHelpers_impl.hpp:69
RsVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &rs)
Constructor.
Definition EquilibrationHelpers_impl.hpp:58
Type that implements "vaporized oil-gas ratio" tabulated as a function of depth policy.
Definition EquilibrationHelpers.hpp:319
Scalar operator()(const Scalar depth, const Scalar press, const Scalar temp, const Scalar satOil=0.0) const override
Function call.
Definition EquilibrationHelpers_impl.hpp:181
RvVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &rv)
Constructor.
Definition EquilibrationHelpers_impl.hpp:170
Type that implements "vaporized water-gas ratio" tabulated as a function of depth policy.
Definition EquilibrationHelpers.hpp:370
Scalar operator()(const Scalar depth, const Scalar press, const Scalar temp, const Scalar satWat=0.0) const override
Function call.
Definition EquilibrationHelpers_impl.hpp:224
RvwVD(const int pvtRegionIdx, const std::vector< Scalar > &depth, const std::vector< Scalar > &rvw)
Constructor.
Definition EquilibrationHelpers_impl.hpp:213
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)
Compute saturation of some phase corresponding to a given capillary pressure.
Definition EquilibrationHelpers_impl.hpp:586
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)
Compute saturation of some phase corresponding to a given capillary pressure, where the capillary pre...
Definition EquilibrationHelpers_impl.hpp:619
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
Functor for inverting a sum of capillary pressure functions.
Definition EquilibrationHelpers.hpp:773
Functor for inverting capillary pressure function.
Definition EquilibrationHelpers.hpp:732