diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index e97d5077..9b062664 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/opm/core/simulator/EquilibrationHelpers.hpp @@ -745,6 +745,7 @@ namespace Opm const PcEq f(props, phase, cell, target_pc); const double f0 = f(s0); const double f1 = f(s1); + if (f0 <= 0.0) { return s0; } else if (f1 > 0.0) { @@ -833,6 +834,46 @@ namespace Opm } } + /// Compute saturation from depth. Used for constant capillary pressure function + inline double satFromDepth(const BlackoilPropertiesInterface& props, + const double cellDepth, + const double contactDepth, + const int phase, + const int cell, + const bool increasing = false) + { + // Find minimum and maximum saturations. + double sminarr[BlackoilPhases::MaxNumPhases]; + double smaxarr[BlackoilPhases::MaxNumPhases]; + props.satRange(1, &cell, sminarr, smaxarr); + const double s0 = increasing ? smaxarr[phase] : sminarr[phase]; + const double s1 = increasing ? sminarr[phase] : smaxarr[phase]; + + if (cellDepth < contactDepth){ + return s0; + } else { + return s1; + } + + } + + /// Return true if capillary pressure function is constant + bool isConstPc(const BlackoilPropertiesInterface& props, + const int phase, + const int cell) + { + // Find minimum and maximum saturations. + double sminarr[BlackoilPhases::MaxNumPhases]; + double smaxarr[BlackoilPhases::MaxNumPhases]; + props.satRange(1, &cell, sminarr, smaxarr); + + // Create the equation f(s) = pc(s); + const PcEq f(props, phase, cell, 0); + const double f0 = f(sminarr[phase]); + const double f1 = f(smaxarr[phase]); + return std::abs(f0 - f1 < std::numeric_limits::epsilon()); + } + } // namespace Equil } // namespace Opm