Make compressible pressure solver compatible with well management.

This commit is contained in:
Xavier Raynaud 2012-05-10 12:38:29 +02:00
parent 2df84c74a8
commit 1e7101a91b
3 changed files with 28 additions and 9 deletions

View File

@ -324,7 +324,11 @@ main(int argc, char** argv)
// Gravity. // Gravity.
gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity; gravity[2] = deck.hasField("NOGRAV") ? 0.0 : Opm::unit::gravity;
// Init state variables (saturation and pressure). // Init state variables (saturation and pressure).
if (param.has("init_saturation")) {
initStateTwophaseBasic(*grid->c_grid(), *props, param, gravity[2], state);
} else {
initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); initStateTwophaseFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
}
} else { } else {
// Grid init. // Grid init.
const int nx = param.getDefault("nx", 100); const int nx = param.getDefault("nx", 100);
@ -518,10 +522,11 @@ main(int argc, char** argv)
std::vector<double> well_perfrates; std::vector<double> well_perfrates;
std::vector<double> fractional_flows; std::vector<double> fractional_flows;
std::vector<double> well_resflows_phase; std::vector<double> well_resflows_phase;
int num_wells = 0;
if (wells->c_wells()) { if (wells->c_wells()) {
const int nw = wells->c_wells()->number_of_wells; num_wells = wells->c_wells()->number_of_wells;
well_bhp.resize(nw, 0.0); well_bhp.resize(num_wells, 0.0);
well_perfrates.resize(wells->c_wells()->well_connpos[nw], 0.0); well_perfrates.resize(wells->c_wells()->well_connpos[num_wells], 0.0);
well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0);
wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates); wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates);
} }
@ -557,15 +562,17 @@ main(int argc, char** argv)
std::vector<double> initial_pressure = state.pressure(); std::vector<double> initial_pressure = state.pressure();
std::vector<double> initial_porevolume(num_cells); std::vector<double> initial_porevolume(num_cells);
computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume); computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume);
std::vector<double> pressure_increment(num_cells); std::vector<double> pressure_increment(num_cells + num_wells);
std::vector<double> prev_pressure; std::vector<double> prev_pressure(num_cells + num_wells);
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
for (int cell = 0; cell < num_cells; ++cell) { for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp->rockComp(state.pressure()[cell]); rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
} }
computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
prev_pressure = state.pressure(); std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin());
std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells);
// prev_pressure = state.pressure();
// compute pressure increment // compute pressure increment
psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc,
@ -577,6 +584,10 @@ main(int argc, char** argv)
state.pressure()[cell] += pressure_increment[cell]; state.pressure()[cell] += pressure_increment[cell];
max_change = std::max(max_change, std::fabs(pressure_increment[cell])); max_change = std::max(max_change, std::fabs(pressure_increment[cell]));
} }
for (int well = 0; well < num_wells; ++well) {
well_bhp[well] += pressure_increment[num_cells + well];
max_change = std::max(max_change, std::fabs(pressure_increment[num_cells + well]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) { if (max_change < nl_pressure_tolerance) {

View File

@ -336,7 +336,10 @@ namespace Opm
soln.well_press = &well_bhp[0]; soln.well_press = &well_bhp[0];
} }
memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x)); // memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x));
ASSERT(int(pressure.size()) == grid_.number_of_cells);
std::copy(pressure.begin(), pressure.end(), h_->x);
std::copy(well_bhp.begin(), well_bhp.end(), h_->x + grid_.number_of_cells);
ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln); ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln);
} }

View File

@ -760,7 +760,7 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G ,
struct ifs_tpfa_data *h ) struct ifs_tpfa_data *h )
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ {
int c, system_singular, ok; int c, w, system_singular, ok;
size_t j; size_t j;
double *v; double *v;
@ -780,6 +780,11 @@ ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G ,
h->A->sa[j] += porevol[c] * rock_comp[c] / dt; h->A->sa[j] += porevol[c] * rock_comp[c] / dt;
h->b[c] += -(porevol[c] - initial_porevolume[c])/dt - v[c]; h->b[c] += -(porevol[c] - initial_porevolume[c])/dt - v[c];
} }
if (F->W != NULL) {
for (w = 0; w < F->W->number_of_wells; ++w) {
h->b[c] += -v[c];
}
}
} }
return ok; return ok;