Use the average well block pressure when evaluating the properties

The average well block pressure is used instead of the well cell
pressure when the well properties are evaluated.
Temperature, rs, rv, phase conditions are still well cells values.

Perforation pressures are stored in the well state
This commit is contained in:
Tor Harald Sandve 2015-02-19 09:50:41 +01:00
parent 46cfaa1432
commit 70a6a2ebf8
2 changed files with 54 additions and 23 deletions

View File

@ -630,10 +630,27 @@ namespace detail {
// Note that some of the complexity of this part is due to the function // Note that some of the complexity of this part is due to the function
// taking std::vector<double> arguments, and not Eigen objects. // taking std::vector<double> arguments, and not Eigen objects.
const int nperf = wells().well_connpos[wells().number_of_wells]; const int nperf = wells().well_connpos[wells().number_of_wells];
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf); const std::vector<int> well_cells(wells().well_cells, wells().well_cells + nperf);
// Compute the average pressure in each well block
const V perf_press = Eigen::Map<const V>(xw.perfPress().data(), nperf);
V avg_press = perf_press*0;
for (int w = 0; w < nw; ++w) {
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) {
const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1];
const double p_avg = (perf_press[perf] + p_above)/2;
avg_press[perf] = p_avg;
}
}
// Use cell values for the temperature as the wells don't knows its temperature yet.
const V perf_temp = subset(state.temperature.value(), well_cells);
// Compute b, rsmax, rvmax values for perforations. // Compute b, rsmax, rvmax values for perforations.
const std::vector<ADB> pressures = computePressures(state); // Evaluate the properties using average well block pressures
const ADB perf_temp = subset(state.temperature, well_cells); // and cell values for rs, rv, phase condition and temperature.
std::vector<PhasePresence> perf_cond(nperf); std::vector<PhasePresence> perf_cond(nperf);
const std::vector<PhasePresence>& pc = phaseCondition(); const std::vector<PhasePresence>& pc = phaseCondition();
for (int perf = 0; perf < nperf; ++perf) { for (int perf = 0; perf < nperf; ++perf) {
@ -641,30 +658,27 @@ namespace detail {
} }
const PhaseUsage& pu = fluid_.phaseUsage(); const PhaseUsage& pu = fluid_.phaseUsage();
DataBlock b(nperf, pu.num_phases); DataBlock b(nperf, pu.num_phases);
std::vector<double> rssat_perf(nperf, 0.0); std::vector<double> rsmax_perf(nperf, 0.0);
std::vector<double> rvsat_perf(nperf, 0.0); std::vector<double> rvmax_perf(nperf, 0.0);
if (pu.phase_used[BlackoilPhases::Aqua]) { if (pu.phase_used[BlackoilPhases::Aqua]) {
const ADB perf_press = subset(pressures[ Water ], well_cells); const V bw = fluid_.bWat(avg_press, perf_temp, well_cells);
const ADB bw = fluid_.bWat(perf_press, perf_temp, well_cells); b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw;
b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value();
} }
assert(active_[Oil]); assert(active_[Oil]);
const ADB perf_so = subset(state.saturation[pu.phase_pos[Oil]], well_cells); const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells);
if (pu.phase_used[BlackoilPhases::Liquid]) { if (pu.phase_used[BlackoilPhases::Liquid]) {
const ADB perf_rs = subset(state.rs, well_cells); const V perf_rs = subset(state.rs.value(), well_cells);
const ADB perf_press = subset(pressures[ Oil ], well_cells); const V bo = fluid_.bOil(avg_press, perf_temp, perf_rs, perf_cond, well_cells);
const ADB bo = fluid_.bOil(perf_press, perf_temp, perf_rs, perf_cond, well_cells); b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo;
b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value(); const V rssat = fluidRsSat(avg_press, perf_so, well_cells);
const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells); rsmax_perf.assign(rssat.data(), rssat.data() + nperf);
rssat_perf.assign(rssat.data(), rssat.data() + nperf);
} }
if (pu.phase_used[BlackoilPhases::Vapour]) { if (pu.phase_used[BlackoilPhases::Vapour]) {
const ADB perf_rv = subset(state.rv, well_cells); const V perf_rv = subset(state.rv.value(), well_cells);
const ADB perf_press = subset(pressures[ Gas ], well_cells); const V bg = fluid_.bGas(avg_press, perf_temp, perf_rv, perf_cond, well_cells);
const ADB bg = fluid_.bGas(perf_press, perf_temp, perf_rv, perf_cond, well_cells); b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value(); const V rvsat = fluidRvSat(avg_press, perf_so, well_cells);
const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells); rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf);
rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf);
} }
// b is row major, so can just copy data. // b is row major, so can just copy data.
std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases); std::vector<double> b_perf(b.data(), b.data() + nperf * pu.num_phases);
@ -689,7 +703,7 @@ namespace detail {
// 2. Compute pressure deltas, and store the results. // 2. Compute pressure deltas, and store the results.
std::vector<double> cdp = WellDensitySegmented std::vector<double> cdp = WellDensitySegmented
::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(), ::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(),
b_perf, rssat_perf, rvsat_perf, perf_depth, b_perf, rsmax_perf, rvmax_perf, perf_depth,
surf_dens, grav); surf_dens, grav);
well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf); well_perforation_pressure_diffs_ = Eigen::Map<const V>(cdp.data(), nperf);
} }
@ -811,8 +825,13 @@ namespace detail {
// DUMPVAL(state.bhp); // DUMPVAL(state.bhp);
// DUMPVAL(ADB::constant(cdp)); // DUMPVAL(ADB::constant(cdp));
// Perforation pressure
const ADB perfpressure = (wops_.w2p * state.bhp) + cdp;
std::vector<double> perfpressure_d(perfpressure.value().data(), perfpressure.value().data() + nperf);
xw.perfPress() = perfpressure_d;
// Pressure drawdown (also used to determine direction of flow) // Pressure drawdown (also used to determine direction of flow)
const ADB drawdown = p_perfcell - (wops_.w2p * state.bhp + cdp); const ADB drawdown = p_perfcell - perfpressure;
// current injecting connections // current injecting connections
auto connInjInx = drawdown.value() < 0; auto connInjInx = drawdown.value() < 0;

View File

@ -48,6 +48,7 @@ namespace Opm
using BaseType :: wellRates; using BaseType :: wellRates;
using BaseType :: bhp; using BaseType :: bhp;
using BaseType :: perfPress;
/// Allocate and initialize if wells is non-null. Also tries /// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates() /// to give useful initial values to the bhp(), wellRates()
@ -96,6 +97,7 @@ namespace Opm
for (int p = 0; p < np; ++p) { for (int p = 0; p < np; ++p) {
perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well); perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
} }
perfPress()[perf] = state.pressure()[wells->well_cells[perf]];
} }
} }
} }
@ -132,13 +134,14 @@ namespace Opm
} }
// perfPhaseRates // perfPhaseRates
int oldPerf = (*it).second[ 1 ] * np; int oldPerf_idx = (*it).second[ 1 ];
const int num_perf_old_well = (*it).second[ 2 ]; const int num_perf_old_well = (*it).second[ 2 ];
const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex]; const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
// copy perforation rates when the number of perforations is equal, // copy perforation rates when the number of perforations is equal,
// otherwise initialize perfphaserates to well rates divided by the number of perforations. // otherwise initialize perfphaserates to well rates divided by the number of perforations.
if( num_perf_old_well == num_perf_this_well ) if( num_perf_old_well == num_perf_this_well )
{ {
int oldPerf = oldPerf_idx *np;
for (int perf = wells->well_connpos[ newIndex ]*np; for (int perf = wells->well_connpos[ newIndex ]*np;
perf < wells->well_connpos[ newIndex + 1]*np; ++perf, ++oldPerf ) perf < wells->well_connpos[ newIndex + 1]*np; ++perf, ++oldPerf )
{ {
@ -151,6 +154,15 @@ namespace Opm
} }
} }
} }
// perfPressures
if( num_perf_old_well == num_perf_this_well )
{
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{
perfPress()[ perf ] = prevState.perfPress()[ oldPerf_idx ];
}
}
// currentControls // currentControls
const int old_control_index = prevState.currentControls()[ oldIndex ]; const int old_control_index = prevState.currentControls()[ oldIndex ];