OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-18 21:43:27 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1027 from totto82/fix_2p_equil_init

Browse Source 
Fix two phase equil initialization
This commit is contained in:
Atgeirr Flø Rasmussen 2016-06-24 14:42:48 +02:00 committed by GitHub
commit 24fa367fda
2 changed files with 62 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
opm/core/props/satfunc/SaturationPropsFromDeck.cpp
View File

@ -136,9 +136,11 @@ namespace Opm
double* pc, double* pc,
double* dpcds) const double* dpcds) const
{ {
assert(cells != 0); assert(cells != 0);
assert(phaseUsage_.phase_used[BlackoilPhases::Liquid]);
const int np = numPhases(); const int np = numPhases();
if (dpcds) { if (dpcds) {
ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_); ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_);
typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation; typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation;
@ -151,12 +153,26 @@ namespace Opm
MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState); MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
// copy the values calculated using opm-material to the target arrays // copy the values calculated using opm-material to the target arrays
for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) { for (int canonicalPhaseIdx = 0; canonicalPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalPhaseIdx) {
double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0; // skip unused phases
pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].value; if ( ! phaseUsage_.phase_used[canonicalPhaseIdx]) {
continue;
}
const int pcPhaseIdx = phaseUsage_.phase_pos[canonicalPhaseIdx];
for (int satPhaseIdx = 0; satPhaseIdx < np; ++satPhaseIdx) const double sign = (canonicalPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].derivatives[satPhaseIdx]; // in opm-material the wetting phase is the reference phase
// for two-phase problems i.e water for oil-water system,
// but for flow it is always oil. Add oil (liquid) capillary pressure value
// to shift the reference phase to oil
pc[np*i + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid].value + sign * capillaryPressures[canonicalPhaseIdx].value;
for (int canonicalSatPhaseIdx = 0; canonicalSatPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalSatPhaseIdx) {
if ( ! phaseUsage_.phase_used[canonicalSatPhaseIdx])
continue;
const int satPhaseIdx = phaseUsage_.phase_pos[canonicalSatPhaseIdx];
dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid].derivatives[canonicalSatPhaseIdx] + sign * capillaryPressures[canonicalPhaseIdx].derivatives[canonicalSatPhaseIdx];
}
} }
} }
} else { } else {
@ -170,9 +186,18 @@ namespace Opm
MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState); MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
// copy the values calculated using opm-material to the target arrays // copy the values calculated using opm-material to the target arrays
for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) { for (int canonicalPhaseIdx = 0; canonicalPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalPhaseIdx) {
double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0; // skip unused phases
pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx]; if ( ! phaseUsage_.phase_used[canonicalPhaseIdx])
continue;
const int pcPhaseIdx = phaseUsage_.phase_pos[canonicalPhaseIdx];
double sign = (canonicalPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
// in opm-material the wetting phase is the reference phase
// for two-phase problems i.e water for oil-water system,
// but for flow it is always oil. Add oil (liquid) capillary pressure value
// to shift the reference phase to oil
pc[np*i + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid] + sign * capillaryPressures[canonicalPhaseIdx];
} }
} }
} }

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

@ -767,28 +767,34 @@ namespace Opm
double sat[BlackoilPhases::MaxNumPhases]; double sat[BlackoilPhases::MaxNumPhases];
double threshold_sat = 1.0e-6; double threshold_sat = 1.0e-6;
sat[waterpos] = smax[waterpos]; sat[oilpos] = 1.0;
sat[gaspos] = smax[gaspos]; if (water) {
sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos]; sat[waterpos] = smax[waterpos];
if (sw > smax[waterpos]-threshold_sat ) { sat[oilpos] -= sat[waterpos];
sat[waterpos] = smax[waterpos]; }
props.capPress(1, sat, &cell, pc, 0); if (gas) {
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos]; sat[gaspos] = smax[gaspos];
} else if (sg > smax[gaspos]-threshold_sat) { sat[oilpos] -= sat[gaspos];
sat[gaspos] = smax[gaspos]; }
props.capPress(1, sat, &cell, pc, 0); if (water && sw > smax[waterpos]-threshold_sat ) {
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos]; sat[waterpos] = smax[waterpos];
} props.capPress(1, sat, &cell, pc, 0);
if (sg < smin[gaspos]+threshold_sat) { phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
sat[gaspos] = smin[gaspos]; } else if (gas && sg > smax[gaspos]-threshold_sat) {
props.capPress(1, sat, &cell, pc, 0); sat[gaspos] = smax[gaspos];
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos]; props.capPress(1, sat, &cell, pc, 0);
} phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
if (sw < smin[waterpos]+threshold_sat) { }
sat[waterpos] = smin[waterpos]; if (gas && sg < smin[gaspos]+threshold_sat) {
props.capPress(1, sat, &cell, pc, 0); sat[gaspos] = smin[gaspos];
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos]; props.capPress(1, sat, &cell, pc, 0);
} phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
}
if (water && sw < smin[waterpos]+threshold_sat) {
sat[waterpos] = smin[waterpos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
}
} }
return phase_saturations; return phase_saturations;
} }