OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

updated polymer_reorder (profiling branch).

Browse Source
This commit is contained in:
Xavier Raynaud 2012-06-13 15:44:34 +02:00
parent a5de5341c5
commit ce6e21c0d6
1 changed files with 42 additions and 80 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

122
examples/polymer_reorder.cpp
View File

@ -28,7 +28,7 @@
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/GridManager.hpp> #include <opm/core/GridManager.hpp>
#include <opm/core/newwells.h> #include <opm/core/newwells.h>
#include <opm/core/WellsManager.hpp> #include <opm/core/wells/WellsManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/initState.hpp> #include <opm/core/utility/initState.hpp>
#include <opm/core/utility/SimulatorTimer.hpp> #include <opm/core/utility/SimulatorTimer.hpp>
@ -53,9 +53,10 @@
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp> #include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp> #include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/ColumnExtract.hpp> #include <opm/core/utility/ColumnExtract.hpp>
#include <opm/polymer/PolymerState.hpp> #include <opm/polymer/PolymerState.hpp>
#include <opm/core/simulator/WellState.hpp>
#include <opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp> #include <opm/polymer/SinglePointUpwindTwoPhasePolymer.hpp>
#include <opm/polymer/GravityColumnSolverPolymer.hpp> #include <opm/polymer/GravityColumnSolverPolymer.hpp>
#include <opm/polymer/TransportModelPolymer.hpp> #include <opm/polymer/TransportModelPolymer.hpp>
@ -570,20 +571,19 @@ main(int argc, char** argv)
// Check that rock compressibility is not used with solvers that do not handle it. // Check that rock compressibility is not used with solvers that do not handle it.
int nl_pressure_maxiter = 0; int nl_pressure_maxiter = 0;
double nl_pressure_tolerance = 0.0; double nl_pressure_residual_tolerance = 0.0;
double nl_pressure_change_tolerance = 0.0;
if (rock_comp->isActive()) { if (rock_comp->isActive()) {
if (!use_reorder) { if (!use_reorder) {
THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet."); THROW("Cannot run implicit (non-reordering) transport solver with rock compressibility yet.");
} }
nl_pressure_residual_tolerance = param.getDefault("nl_pressure_residual_tolerance", 0.0);
nl_pressure_change_tolerance = param.getDefault("nl_pressure_change_tolerance", 1.0); // In Pascal.
nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10); nl_pressure_maxiter = param.getDefault("nl_pressure_maxiter", 10);
nl_pressure_tolerance = param.getDefault("nl_pressure_tolerance", 1.0); // in Pascal
} }
// Source-related variables init. // Source-related variables init.
int num_cells = grid->c_grid()->number_of_cells; int num_cells = grid->c_grid()->number_of_cells;
std::vector<double> totmob;
std::vector<double> omega; // Will remain empty if no gravity.
std::vector<double> rc; // Will remain empty if no rock compressibility.
// Extra rock init. // Extra rock init.
std::vector<double> porevol; std::vector<double> porevol;
@ -594,12 +594,8 @@ main(int argc, char** argv)
} }
double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0); double tot_porevol_init = std::accumulate(porevol.begin(), porevol.end(), 0.0);
// We need a separate reorder_sat, because the reorder
// code expects a scalar sw, not both sw and so.
std::vector<double> reorder_sat(num_cells);
std::vector<double> src(num_cells, 0.0);
// Initialising src // Initialising src
std::vector<double> src(num_cells, 0.0);
if (wells->c_wells()) { if (wells->c_wells()) {
// Do nothing, wells will be the driving force, not source terms. // Do nothing, wells will be the driving force, not source terms.
// Opm::wellsToSrc(*wells->c_wells(), num_cells, src); // Opm::wellsToSrc(*wells->c_wells(), num_cells, src);
@ -625,7 +621,10 @@ main(int argc, char** argv)
Opm::LinearSolverFactory linsolver(param); Opm::LinearSolverFactory linsolver(param);
// Pressure solver. // Pressure solver.
const double *grav = use_gravity ? &gravity[0] : 0; const double *grav = use_gravity ? &gravity[0] : 0;
Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), grav, linsolver, wells->c_wells()); Opm::IncompTpfa psolver(*grid->c_grid(), *props, rock_comp.get(), linsolver,
nl_pressure_residual_tolerance, nl_pressure_change_tolerance,
nl_pressure_maxiter,
grav, wells->c_wells(), src, bcs.c_bcs());
// Reordering solver. // Reordering solver.
const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9); const double nl_tolerance = param.getDefault("nl_tolerance", 1e-9);
const int nl_maxiter = param.getDefault("nl_maxiter", 30); const int nl_maxiter = param.getDefault("nl_maxiter", 30);
@ -639,7 +638,7 @@ main(int argc, char** argv)
THROW("Unknown method: " << method_string); THROW("Unknown method: " << method_string);
} }
Opm::TransportModelPolymer reorder_model(*grid->c_grid(), props->porosity(), &porevol[0], *props, polyprop, Opm::TransportModelPolymer reorder_model(*grid->c_grid(), *props, polyprop,
method, nl_tolerance, nl_maxiter); method, nl_tolerance, nl_maxiter);
if (use_gauss_seidel_gravity) { if (use_gauss_seidel_gravity) {
@ -717,17 +716,15 @@ main(int argc, char** argv)
Opm::Watercut watercut; Opm::Watercut watercut;
watercut.push(0.0, 0.0, 0.0); watercut.push(0.0, 0.0, 0.0);
Opm::WellReport wellreport; Opm::WellReport wellreport;
std::vector<double> well_bhp; Opm::WellState well_state;
std::vector<double> well_perfrates; well_state.init(wells->c_wells(), state);
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; int num_wells = 0;
if (wells->c_wells()) { if (wells->c_wells()) {
num_wells = wells->c_wells()->number_of_wells; num_wells = wells->c_wells()->number_of_wells;
well_bhp.resize(num_wells, 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_state.bhp(), well_state.perfRates());
} }
for (; !simtimer.done(); ++simtimer) { for (; !simtimer.done(); ++simtimer) {
// Report timestep and (optionally) write state to disk. // Report timestep and (optionally) write state to disk.
@ -737,17 +734,6 @@ main(int argc, char** argv)
} }
// Solve pressure. // Solve pressure.
if (use_gravity) {
computeTotalMobilityOmega(*props, polyprop, allcells, state.saturation(), state.concentration(), state.maxconcentration(),
totmob, omega);
} else {
computeTotalMobility(*props, polyprop, allcells, state.saturation(), state.concentration(), state.maxconcentration(),
totmob);
}
std::vector<double> wdp;
if (wells->c_wells()) {
Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), gravity[2], true, wdp);
}
if (check_well_controls) { if (check_well_controls) {
computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows); computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows);
} }
@ -758,48 +744,27 @@ main(int argc, char** argv)
int well_control_iteration = 0; int well_control_iteration = 0;
do { do {
pressure_timer.start(); pressure_timer.start();
if (rock_comp->isActive()) { std::vector<double> initial_pressure = state.pressure();
rc.resize(num_cells); psolver.solve(simtimer.currentStepLength(), state, well_state);
std::vector<double> initial_pressure = state.pressure(); if (!rock_comp->isActive()) {
std::vector<double> initial_porevolume(num_cells); // Compute average pressures of previous and last
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, initial_pressure, initial_porevolume); // step, and total volume.
std::vector<double> pressure_increment(num_cells + num_wells); double av_prev_press = 0.;
std::vector<double> prev_pressure(num_cells + num_wells); double av_press = 0.;
for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { double tot_vol = 0.;
for (int cell = 0; cell < num_cells; ++cell) {
for (int cell = 0; cell < num_cells; ++cell) { av_prev_press += initial_pressure[cell]*grid->c_grid()->cell_volumes[cell];
rc[cell] = rock_comp->rockComp(state.pressure()[cell]); av_press += state.pressure()[cell]*grid->c_grid()->cell_volumes[cell];
} tot_vol += grid->c_grid()->cell_volumes[cell];
computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); }
std::copy(state.pressure().begin(), state.pressure().end(), prev_pressure.begin()); // Renormalization constant
std::copy(well_bhp.begin(), well_bhp.end(), prev_pressure.begin() + num_cells); const double ren_const = (av_prev_press - av_press)/tot_vol;
// prev_pressure = state.pressure(); for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += ren_const;
// compute pressure increment }
psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, for (int well = 0; well < num_wells; ++well) {
prev_pressure, initial_porevolume, simtimer.currentStepLength(), well_state.bhp()[well] += ren_const;
pressure_increment);
double max_change = 0.0;
for (int cell = 0; cell < num_cells; ++cell) {
state.pressure()[cell] += 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;
if (max_change < nl_pressure_tolerance) {
break;
}
} }
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);
} }
pressure_timer.stop(); pressure_timer.stop();
double pt = pressure_timer.secsSinceStart(); double pt = pressure_timer.secsSinceStart();
@ -810,11 +775,11 @@ main(int argc, char** argv)
if (check_well_controls) { if (check_well_controls) {
Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), Opm::computePhaseFlowRatesPerWell(*wells->c_wells(),
fractional_flows, fractional_flows,
well_perfrates, well_state.perfRates(),
well_resflows_phase); well_resflows_phase);
std::cout << "Checking well conditions." << std::endl; std::cout << "Checking well conditions." << std::endl;
// For testing we set surface := reservoir // For testing we set surface := reservoir
well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase); well_control_passed = wells->conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
++well_control_iteration; ++well_control_iteration;
if (!well_control_passed && well_control_iteration > max_well_control_iterations) { if (!well_control_passed && well_control_iteration > max_well_control_iterations) {
THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries.");
@ -829,7 +794,7 @@ main(int argc, char** argv)
// Process transport sources (to include bdy terms and well flows). // Process transport sources (to include bdy terms and well flows).
Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0,
wells->c_wells(), well_perfrates, reorder_src); wells->c_wells(), well_state.perfRates(), reorder_src);
// Find inflow rate. // Find inflow rate.
@ -851,19 +816,16 @@ main(int argc, char** argv)
} }
for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) { for (int tr_substep = 0; tr_substep < num_transport_substeps; ++tr_substep) {
if (use_reorder) { if (use_reorder) {
Opm::toWaterSat(state.saturation(), reorder_sat);
reorder_model.solve(&state.faceflux()[0], &porevol[0], &reorder_src[0], stepsize, inflow_c, reorder_model.solve(&state.faceflux()[0], &porevol[0], &reorder_src[0], stepsize, inflow_c,
&reorder_sat[0], &state.concentration()[0], &state.maxconcentration()[0]); state.saturation(), state.concentration(), state.maxconcentration());
Opm::toBothSat(reorder_sat, state.saturation());
Opm::computeInjectedProduced(*props, polyprop, state.saturation(), state.concentration(), state.maxconcentration(), Opm::computeInjectedProduced(*props, polyprop, state.saturation(), state.concentration(), state.maxconcentration(),
reorder_src, simtimer.currentStepLength(), inflow_c, reorder_src, simtimer.currentStepLength(), inflow_c,
injected, produced, polyinj, polyprod); injected, produced, polyinj, polyprod);
if (use_segregation_split) { if (use_segregation_split) {
if (use_column_solver) { if (use_column_solver) {
if (use_gauss_seidel_gravity) { if (use_gauss_seidel_gravity) {
reorder_model.solveGravity(columns, &porevol[0], stepsize, reorder_sat, reorder_model.solveGravity(columns, &porevol[0], stepsize, state.saturation(),
state.concentration(), state.maxconcentration()); state.concentration(), state.maxconcentration());
Opm::toBothSat(reorder_sat, state.saturation());
} else { } else {
colsolver.solve(columns, stepsize, state.saturation(), state.concentration(), colsolver.solve(columns, stepsize, state.saturation(), state.concentration(),
state.maxconcentration()); state.maxconcentration());
@ -936,7 +898,7 @@ main(int argc, char** argv)
if (wells->c_wells()) { if (wells->c_wells()) {
wellreport.push(*props, *wells->c_wells(), state.saturation(), wellreport.push(*props, *wells->c_wells(), state.saturation(),
simtimer.currentTime() + simtimer.currentStepLength(), simtimer.currentTime() + simtimer.currentStepLength(),
well_bhp, well_perfrates); well_state.bhp(), well_state.perfRates());
} }
} }
total_timer.stop(); total_timer.stop();