21#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
22#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
24#include <opm/core/utility/linearInterpolation.hpp>
25#include <opm/core/utility/RegionMapping.hpp>
26#include <opm/core/utility/RootFinders.hpp>
28#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
30#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
31#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
32#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
93 typedef Opm::SimpleModularFluidState<double,
110 namespace Miscibility {
136 const double sat = 0.0)
const = 0;
165 const double = 0.0)
const
177 template <
class Flu
idSystem>
188 const std::vector<double>& depth,
189 const std::vector<double>& rs)
190 : pvtRegionIdx_(pvtRegionIdx)
215 const double sat_gas = 0.0)
const
218 return satRs(press, temp);
220 return std::min(satRs(press, temp), linearInterpolationNoExtrapolation(depth_, rs_, depth));
225 const int pvtRegionIdx_;
226 std::vector<double> depth_;
227 std::vector<double> rs_;
229 double satRs(
const double press,
const double temp)
const
231 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
241 template <
class Flu
idSystem>
252 const std::vector<double>& depth,
253 const std::vector<double>& rv)
254 : pvtRegionIdx_(pvtRegionIdx)
279 const double sat_oil = 0.0 )
const
281 if (std::abs(sat_oil) > 1e-16) {
282 return satRv(press, temp);
284 return std::min(satRv(press, temp), linearInterpolationNoExtrapolation(depth_, rv_, depth));
289 const int pvtRegionIdx_;
290 std::vector<double> depth_;
291 std::vector<double> rv_;
293 double satRv(
const double press,
const double temp)
const
295 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
314 template <
class Flu
idSystem>
324 RsSatAtContact(
const int pvtRegionIdx,
const double p_contact,
const double T_contact)
325 : pvtRegionIdx_(pvtRegionIdx)
327 rs_sat_contact_ = satRs(p_contact, T_contact);
349 const double sat_gas = 0.0)
const
352 return satRs(press, temp);
354 return std::min(satRs(press, temp), rs_sat_contact_);
359 const int pvtRegionIdx_;
360 double rs_sat_contact_;
362 double satRs(
const double press,
const double temp)
const
364 return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
383 template <
class Flu
idSystem>
393 RvSatAtContact(
const int pvtRegionIdx,
const double p_contact,
const double T_contact)
394 :pvtRegionIdx_(pvtRegionIdx)
396 rv_sat_contact_ = satRv(p_contact, T_contact);
418 const double sat_oil = 0.0)
const
421 return satRv(press, temp);
423 return std::min(satRv(press, temp), rv_sat_contact_);
428 const int pvtRegionIdx_;
429 double rv_sat_contact_;
431 double satRv(
const double press,
const double temp)
const
433 return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
469 std::shared_ptr<Miscibility::RsFunction> rs,
470 std::shared_ptr<Miscibility::RsFunction> rv,
492 double datum()
const {
return this->rec_.datumDepth(); }
497 double pressure()
const {
return this->rec_.datumDepthPressure(); }
502 double zwoc()
const {
return this->rec_.waterOilContactDepth(); }
509 double pcow_woc()
const {
return this->rec_.waterOilContactCapillaryPressure(); }
514 double zgoc()
const {
return this->rec_.gasOilContactDepth(); }
521 double pcgo_goc()
const {
return this->rec_.gasOilContactCapillaryPressure(); }
541 int pvtIdx()
const {
return this->pvtIdx_; }
546 std::shared_ptr<Miscibility::RsFunction> rs_;
547 std::shared_ptr<Miscibility::RsFunction> rv_;
556 template <
class Flu
idSystem,
class MaterialLaw,
class MaterialLawManager>
559 PcEq(
const MaterialLawManager& materialLawManager,
562 const double target_pc)
563 : materialLawManager_(materialLawManager),
566 target_pc_(target_pc)
572 const auto& matParams = materialLawManager_.materialLawParams(cell_);
574 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
575 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
576 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
577 fluidState.setSaturation(phase_, s);
579 double pc[FluidSystem::numPhases];
580 std::fill(pc, pc + FluidSystem::numPhases, 0.0);
581 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
582 double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
583 double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign * pc[phase_];
584 return pcPhase - target_pc_;
587 const MaterialLawManager& materialLawManager_;
590 const double target_pc_;
593 template <
class Flu
idSystem,
class MaterialLawManager>
594 double minSaturations(
const MaterialLawManager& materialLawManager,
const int phase,
const int cell) {
595 const auto& scaledDrainageInfo =
596 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
600 case FluidSystem::waterPhaseIdx :
602 return scaledDrainageInfo.Swl;
605 case FluidSystem::gasPhaseIdx :
607 return scaledDrainageInfo.Sgl;
610 case FluidSystem::oilPhaseIdx :
612 OPM_THROW(std::runtime_error,
"Min saturation not implemented for oil phase.");
615 default: OPM_THROW(std::runtime_error,
"Unknown phaseIdx .");
620 template <
class Flu
idSystem,
class MaterialLawManager>
621 double maxSaturations(
const MaterialLawManager& materialLawManager,
const int phase,
const int cell) {
622 const auto& scaledDrainageInfo =
623 materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
627 case FluidSystem::waterPhaseIdx :
629 return scaledDrainageInfo.Swu;
632 case FluidSystem::gasPhaseIdx :
634 return scaledDrainageInfo.Sgu;
637 case FluidSystem::oilPhaseIdx :
639 OPM_THROW(std::runtime_error,
"Max saturation not implemented for oil phase.");
642 default: OPM_THROW(std::runtime_error,
"Unknown phaseIdx .");
650 template <
class Flu
idSystem,
class MaterialLaw,
class MaterialLawManager >
651 inline double satFromPc(
const MaterialLawManager& materialLawManager,
654 const double target_pc,
655 const bool increasing =
false)
658 const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
659 const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
663 const double f0 = f(s0);
664 const double f1 = f(s1);
668 }
else if (f1 > 0.0) {
671 const int max_iter = 60;
672 const double tol = 1e-6;
674 typedef RegulaFalsi<ThrowOnError> ScalarSolver;
675 const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used);
684 template <
class Flu
idSystem,
class MaterialLaw,
class MaterialLawManager>
687 PcEqSum(
const MaterialLawManager& materialLawManager,
691 const double target_pc)
692 : materialLawManager_(materialLawManager),
696 target_pc_(target_pc)
701 const auto& matParams = materialLawManager_.materialLawParams(cell_);
703 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
704 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
705 fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
706 fluidState.setSaturation(phase1_, s);
707 fluidState.setSaturation(phase2_, 1.0 - s);
709 double pc[FluidSystem::numPhases];
710 std::fill(pc, pc + FluidSystem::numPhases, 0.0);
712 MaterialLaw::capillaryPressures(pc, matParams, fluidState);
713 double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
714 double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 * pc[phase1_];
715 double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
716 double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 * pc[phase2_];
717 return pc1 + pc2 - target_pc_;
720 const MaterialLawManager& materialLawManager_;
724 const double target_pc_;
733 template <
class Flu
idSystem,
class MaterialLaw,
class MaterialLawManager>
738 const double target_pc)
741 const double smin = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
742 const double smax = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
746 const double f0 = f(smin);
747 const double f1 = f(smax);
750 }
else if (f1 > 0.0) {
753 const int max_iter = 30;
754 const double tol = 1e-6;
756 typedef RegulaFalsi<ThrowOnError> ScalarSolver;
757 const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
763 template <
class Flu
idSystem,
class MaterialLaw,
class MaterialLawManager>
764 inline double satFromDepth(
const MaterialLawManager& materialLawManager,
765 const double cellDepth,
766 const double contactDepth,
769 const bool increasing =
false)
771 const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
772 const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
774 if (cellDepth < contactDepth){
783 template <
class Flu
idSystem,
class MaterialLaw,
class MaterialLawManager>
784 inline bool isConstPc(
const MaterialLawManager& materialLawManager,
790 const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
791 const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
792 return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
Definition: EquilibrationHelpers.hpp:458
EquilReg(const EquilRecord &rec, std::shared_ptr< Miscibility::RsFunction > rs, std::shared_ptr< Miscibility::RsFunction > rv, const int pvtIdx)
Definition: EquilibrationHelpers.hpp:468
double pcgo_goc() const
Definition: EquilibrationHelpers.hpp:521
double datum() const
Definition: EquilibrationHelpers.hpp:492
double zgoc() const
Definition: EquilibrationHelpers.hpp:514
double pcow_woc() const
Definition: EquilibrationHelpers.hpp:509
Miscibility::RsFunction CalcDissolution
Definition: EquilibrationHelpers.hpp:482
double zwoc() const
Definition: EquilibrationHelpers.hpp:502
Miscibility::RsFunction CalcEvaporation
Definition: EquilibrationHelpers.hpp:487
double pressure() const
Definition: EquilibrationHelpers.hpp:497
int pvtIdx() const
Definition: EquilibrationHelpers.hpp:541
const CalcEvaporation & evaporationCalculator() const
Definition: EquilibrationHelpers.hpp:536
const CalcDissolution & dissolutionCalculator() const
Definition: EquilibrationHelpers.hpp:529
Definition: EquilibrationHelpers.hpp:143
double operator()(const double, const double, const double, const double=0.0) const
Definition: EquilibrationHelpers.hpp:162
Definition: EquilibrationHelpers.hpp:116
virtual double operator()(const double depth, const double press, const double temp, const double sat=0.0) const =0
Definition: EquilibrationHelpers.hpp:178
RsVD(const int pvtRegionIdx, const std::vector< double > &depth, const std::vector< double > &rs)
Definition: EquilibrationHelpers.hpp:187
double operator()(const double depth, const double press, const double temp, const double sat_gas=0.0) const
Definition: EquilibrationHelpers.hpp:212
Definition: EquilibrationHelpers.hpp:242
double operator()(const double depth, const double press, const double temp, const double sat_oil=0.0) const
Definition: EquilibrationHelpers.hpp:276
RvVD(const int pvtRegionIdx, const std::vector< double > &depth, const std::vector< double > &rv)
Definition: EquilibrationHelpers.hpp:251
double satFromSumOfPcs(const MaterialLawManager &materialLawManager, const int phase1, const int phase2, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:734
Opm::FluidSystems::BlackOil< double > FluidSystemSimple
Definition: EquilibrationHelpers.hpp:90
Opm::SimpleModularFluidState< double, 3, 3, FluidSystemSimple, false, false, false, false, true, false, false, false > SatOnlyFluidState
Definition: EquilibrationHelpers.hpp:104
bool isConstPc(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition: EquilibrationHelpers.hpp:784
double satFromPc(const MaterialLawManager &materialLawManager, const int phase, const int cell, const double target_pc, const bool increasing=false)
Definition: EquilibrationHelpers.hpp:651
double maxSaturations(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Definition: EquilibrationHelpers.hpp:621
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.hpp:764
double minSaturations(const MaterialLawManager &materialLawManager, const int phase, const int cell)
Definition: EquilibrationHelpers.hpp:594
Definition: AnisotropicEikonal.hpp:44
Definition: EquilibrationHelpers.hpp:558
double operator()(double s) const
Definition: EquilibrationHelpers.hpp:570
PcEq(const MaterialLawManager &materialLawManager, const int phase, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:559
Definition: EquilibrationHelpers.hpp:686
PcEqSum(const MaterialLawManager &materialLawManager, const int phase1, const int phase2, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:687
double operator()(double s) const
Definition: EquilibrationHelpers.hpp:699