From 4f862e759e9b2038dee62b526be2e6bdd256c583 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 11 Aug 2014 11:21:44 +0200 Subject: [PATCH 1/2] Computes saturations based on depths For constant capillar pressure function the saturation is determined by cell depths: Sg_max, Sw_min ----- goc ---- Sg_min, Sw_min ----- woc ---- Sg_min, Sw_max --- opm/core/simulator/EquilibrationHelpers.hpp | 41 +++++++++++++++++++++ 1 file changed, 41 insertions(+) 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 From 5051658a07b15e3b81526b72be928b2220d25d6c Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 11 Aug 2014 11:23:15 +0200 Subject: [PATCH 2/2] Implements initialization for constant capillary pressure functions --- opm/core/simulator/initStateEquil.hpp | 2 +- opm/core/simulator/initStateEquil_impl.hpp | 45 +++++++++++++++------- 2 files changed, 33 insertions(+), 14 deletions(-) diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index b869aa97..bcee7145 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -403,7 +403,7 @@ namespace Opm PVec press = phasePressures(G, eqreg, cells, grav); - const PVec sat = phaseSaturations(eqreg, cells, props, swat_init_, press); + const PVec sat = phaseSaturations(G, eqreg, cells, props, swat_init_, press); const int np = props.numPhases(); for (int p = 0; p < np; ++p) { diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index 611c5cd4..5e68a176 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -578,7 +578,8 @@ namespace Opm template std::vector< std::vector > - phaseSaturations(const Region& reg, + phaseSaturations(const UnstructuredGrid& G, + const Region& reg, const CellRange& cells, BlackoilPropertiesInterface& props, const std::vector swat_init, @@ -598,6 +599,8 @@ namespace Opm double smin[BlackoilPhases::MaxNumPhases] = { 0.0 }; double smax[BlackoilPhases::MaxNumPhases] = { 0.0 }; + const double* const cc = & G.cell_centroids[0]; + const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua]; const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour]; const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid]; @@ -611,23 +614,39 @@ namespace Opm // inverting capillary pressure functions. double sw = 0.0; if (water) { - const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; - if (swat_init.empty()) { // Invert Pc to find sw - sw = satFromPc(props, waterpos, cell, pcov); - phase_saturations[waterpos][local_index] = sw; - } else { // Scale Pc to reflect imposed sw - sw = swat_init[cell]; - props.swatInitScaling(cell, pcov, sw); + if (isConstPc(props,waterpos,cell)){ + const int nd = G.dimensions; + const double cellDepth = cc[nd * cell + nd-1]; + sw = satFromDepth(props,cellDepth,zwoc,waterpos,cell,false); phase_saturations[waterpos][local_index] = sw; } + else{ + const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; + if (swat_init.empty()) { // Invert Pc to find sw + sw = satFromPc(props, waterpos, cell, pcov); + phase_saturations[waterpos][local_index] = sw; + } else { // Scale Pc to reflect imposed sw + sw = swat_init[cell]; + props.swatInitScaling(cell, pcov, sw); + phase_saturations[waterpos][local_index] = sw; + } + } } double sg = 0.0; if (gas) { - // Note that pcog is defined to be (pg - po), not (po - pg). - const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index]; - const double increasing = true; // pcog(sg) expected to be increasing function - sg = satFromPc(props, gaspos, cell, pcog, increasing); - phase_saturations[gaspos][local_index] = sg; + if (isConstPc(props,gaspos,cell)){ + const int nd = G.dimensions; + const double cellDepth = cc[nd * cell + nd-1]; + sg = satFromDepth(props,cellDepth,zgoc,gaspos,cell,true); + phase_saturations[gaspos][local_index] = sg; + } + else{ + // Note that pcog is defined to be (pg - po), not (po - pg). + const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index]; + const double increasing = true; // pcog(sg) expected to be increasing function + sg = satFromPc(props, gaspos, cell, pcog, increasing); + phase_saturations[gaspos][local_index] = sg; + } } if (gas && water && (sg + sw > 1.0)) { // Overlapping gas-oil and oil-water transition