OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #628 from totto82/constCapPressInit2

Browse Source 
Initialization for constant capillary pressure functions
This commit is contained in:
Atgeirr Flø Rasmussen 2014-08-11 13:49:37 +02:00
commit 1d4e672c61
3 changed files with 74 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
opm/core/simulator/EquilibrationHelpers.hpp
View File

@ -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<double>::epsilon());
}
} // namespace Equil
} // namespace Opm

2
opm/core/simulator/initStateEquil.hpp
View File

@ -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) {

45
opm/core/simulator/initStateEquil_impl.hpp
View File

@ -578,7 +578,8 @@ namespace Opm
template <class Region, class CellRange>
std::vector< std::vector<double> >
phaseSaturations(const Region& reg,
phaseSaturations(const UnstructuredGrid& G,
const Region& reg,
const CellRange& cells,
BlackoilPropertiesInterface& props,
const std::vector<double> 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