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

Modified spu_2p to agree with newton compressible pressure solver.

Browse Source
This commit is contained in:
Xavier Raynaud 2012-05-09 14:49:23 +02:00
parent a2e5c4b44e
commit 933d4adf44
1 changed files with 24 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
examples/spu_2p.cpp
View File

@ -539,25 +539,39 @@ main(int argc, char** argv)
if (rock_comp->isActive()) {
rc.resize(num_cells);
std::vector<double> initial_pressure = state.pressure();
std::vector<double> initial_porevolume(num_cells);
computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume);
std::vector<double> pressure_increment(num_cells);
std::vector<double> prev_pressure;
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
prev_pressure = state.pressure();
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
}
state.pressure() = initial_pressure;
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(),
state.pressure(), state.faceflux(), well_bhp, well_perfrates);
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) {
for (int cell = 0; cell < num_cells; ++cell) {
rc[cell] = rock_comp->rockComp(state.pressure()[cell]);
}
computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
prev_pressure = state.pressure();
// compute pressure increment
psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc,
prev_pressure, initial_porevolume, simtimer.currentStepLength(),
pressure_increment);
double max_change = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell]));
state.pressure()[cell] += pressure_increment[cell];
max_change = std::max(max_change, std::fabs(pressure_increment[cell]));
}
std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl;
if (max_change < nl_pressure_tolerance) {
break;
}
}
computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol);
psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);
} else {
psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(),
well_bhp, well_perfrates);