mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
due to the new timer, simulator class should not do time iteration once again.
This commit is contained in:
Liu Ming
parent
5242d6bbf7
commit
c0a61c9655
@ -292,208 +292,206 @@ namespace Opm
|
|||||||
wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
|
wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
|
||||||
state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
|
state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
|
||||||
}
|
}
|
||||||
// for (; !timer.done(); ++timer) {
|
// Report timestep and (optionally) write state to disk.
|
||||||
// Report timestep and (optionally) write state to disk.
|
timer.report(std::cout);
|
||||||
timer.report(std::cout);
|
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
|
||||||
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
|
if (output_vtk_) {
|
||||||
if (output_vtk_) {
|
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
}
|
||||||
|
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
|
}
|
||||||
|
|
||||||
|
initial_pressure = state.pressure();
|
||||||
|
|
||||||
|
// Solve pressure equation.
|
||||||
|
if (check_well_controls_) {
|
||||||
|
computeFractionalFlow(props_, poly_props_, allcells_,
|
||||||
|
state.pressure(), state.surfacevol(), state.saturation(),
|
||||||
|
state.concentration(), state.maxconcentration(),
|
||||||
|
fractional_flows);
|
||||||
|
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
|
||||||
|
}
|
||||||
|
bool well_control_passed = !check_well_controls_;
|
||||||
|
int well_control_iteration = 0;
|
||||||
|
do {
|
||||||
|
// Run solver
|
||||||
|
pressure_timer.start();
|
||||||
|
psolver_.solve(timer.currentStepLength(), state, well_state);
|
||||||
|
|
||||||
|
// Renormalize pressure if both fluids and rock are
|
||||||
|
// incompressible, and there are no pressure
|
||||||
|
// conditions (bcs or wells). It is deemed sufficient
|
||||||
|
// for now to renormalize using geometric volume
|
||||||
|
// instead of pore volume.
|
||||||
|
if (psolver_.singularPressure()) {
|
||||||
|
// Compute average pressures of previous and last
|
||||||
|
// step, and total volume.
|
||||||
|
double av_prev_press = 0.0;
|
||||||
|
double av_press = 0.0;
|
||||||
|
double tot_vol = 0.0;
|
||||||
|
const int num_cells = grid_.number_of_cells;
|
||||||
|
for (int cell = 0; cell < num_cells; ++cell) {
|
||||||
|
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
|
||||||
|
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
|
||||||
|
tot_vol += grid_.cell_volumes[cell];
|
||||||
|
}
|
||||||
|
// Renormalization constant
|
||||||
|
const double ren_const = (av_prev_press - av_press)/tot_vol;
|
||||||
|
for (int cell = 0; cell < num_cells; ++cell) {
|
||||||
|
state.pressure()[cell] += ren_const;
|
||||||
|
}
|
||||||
|
const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
|
||||||
|
for (int well = 0; well < num_wells; ++well) {
|
||||||
|
well_state.bhp()[well] += ren_const;
|
||||||
}
|
}
|
||||||
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
initial_pressure = state.pressure();
|
// Stop timer and report
|
||||||
|
pressure_timer.stop();
|
||||||
|
double pt = pressure_timer.secsSinceStart();
|
||||||
|
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
|
||||||
|
ptime += pt;
|
||||||
|
|
||||||
// Solve pressure equation.
|
// Optionally, check if well controls are satisfied.
|
||||||
if (check_well_controls_) {
|
if (check_well_controls_) {
|
||||||
computeFractionalFlow(props_, poly_props_, allcells_,
|
Opm::computePhaseFlowRatesPerWell(*wells_,
|
||||||
state.pressure(), state.surfacevol(), state.saturation(),
|
well_state.perfRates(),
|
||||||
state.concentration(), state.maxconcentration(),
|
fractional_flows,
|
||||||
fractional_flows);
|
well_resflows_phase);
|
||||||
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
|
std::cout << "Checking well conditions." << std::endl;
|
||||||
}
|
// For testing we set surface := reservoir
|
||||||
bool well_control_passed = !check_well_controls_;
|
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
|
||||||
int well_control_iteration = 0;
|
++well_control_iteration;
|
||||||
do {
|
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
|
||||||
// Run solver
|
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
|
||||||
pressure_timer.start();
|
|
||||||
psolver_.solve(timer.currentStepLength(), state, well_state);
|
|
||||||
|
|
||||||
// Renormalize pressure if both fluids and rock are
|
|
||||||
// incompressible, and there are no pressure
|
|
||||||
// conditions (bcs or wells). It is deemed sufficient
|
|
||||||
// for now to renormalize using geometric volume
|
|
||||||
// instead of pore volume.
|
|
||||||
if (psolver_.singularPressure()) {
|
|
||||||
// Compute average pressures of previous and last
|
|
||||||
// step, and total volume.
|
|
||||||
double av_prev_press = 0.0;
|
|
||||||
double av_press = 0.0;
|
|
||||||
double tot_vol = 0.0;
|
|
||||||
const int num_cells = grid_.number_of_cells;
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
|
|
||||||
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
|
|
||||||
tot_vol += grid_.cell_volumes[cell];
|
|
||||||
}
|
|
||||||
// Renormalization constant
|
|
||||||
const double ren_const = (av_prev_press - av_press)/tot_vol;
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
state.pressure()[cell] += ren_const;
|
|
||||||
}
|
|
||||||
const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
|
|
||||||
for (int well = 0; well < num_wells; ++well) {
|
|
||||||
well_state.bhp()[well] += ren_const;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
if (!well_control_passed) {
|
||||||
// Stop timer and report
|
std::cout << "Well controls not passed, solving again." << std::endl;
|
||||||
pressure_timer.stop();
|
} else {
|
||||||
double pt = pressure_timer.secsSinceStart();
|
std::cout << "Well conditions met." << std::endl;
|
||||||
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
|
|
||||||
ptime += pt;
|
|
||||||
|
|
||||||
// Optionally, check if well controls are satisfied.
|
|
||||||
if (check_well_controls_) {
|
|
||||||
Opm::computePhaseFlowRatesPerWell(*wells_,
|
|
||||||
well_state.perfRates(),
|
|
||||||
fractional_flows,
|
|
||||||
well_resflows_phase);
|
|
||||||
std::cout << "Checking well conditions." << std::endl;
|
|
||||||
// For testing we set surface := reservoir
|
|
||||||
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
|
|
||||||
++well_control_iteration;
|
|
||||||
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
|
|
||||||
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
|
|
||||||
}
|
|
||||||
if (!well_control_passed) {
|
|
||||||
std::cout << "Well controls not passed, solving again." << std::endl;
|
|
||||||
} else {
|
|
||||||
std::cout << "Well conditions met." << std::endl;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
} while (!well_control_passed);
|
|
||||||
|
|
||||||
// Update pore volumes if rock is compressible.
|
|
||||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
|
||||||
initial_porevol = porevol;
|
|
||||||
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Process transport sources (to include bdy terms and well flows).
|
|
||||||
Opm::computeTransportSource(props_, wells_, well_state, transport_src);
|
|
||||||
|
|
||||||
// Find inflow rate.
|
|
||||||
const double current_time = timer.simulationTimeElapsed();
|
|
||||||
double stepsize = timer.currentStepLength();
|
|
||||||
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
|
|
||||||
|
|
||||||
|
|
||||||
// Solve transport.
|
|
||||||
transport_timer.start();
|
|
||||||
if (num_transport_substeps_ != 1) {
|
|
||||||
stepsize /= double(num_transport_substeps_);
|
|
||||||
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
|
|
||||||
}
|
|
||||||
double injected[2] = { 0.0 };
|
|
||||||
double produced[2] = { 0.0 };
|
|
||||||
double polyinj = 0.0;
|
|
||||||
double polyprod = 0.0;
|
|
||||||
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
|
|
||||||
tsolver_.solve(&state.faceflux()[0], initial_pressure,
|
|
||||||
state.pressure(), &initial_porevol[0], &porevol[0],
|
|
||||||
&transport_src[0], &polymer_inflow_c[0], stepsize,
|
|
||||||
state.saturation(), state.surfacevol(),
|
|
||||||
state.concentration(), state.maxconcentration());
|
|
||||||
double substep_injected[2] = { 0.0 };
|
|
||||||
double substep_produced[2] = { 0.0 };
|
|
||||||
double substep_polyinj = 0.0;
|
|
||||||
double substep_polyprod = 0.0;
|
|
||||||
Opm::computeInjectedProduced(props_, poly_props_,
|
|
||||||
state,
|
|
||||||
transport_src, polymer_inflow_c, stepsize,
|
|
||||||
substep_injected, substep_produced,
|
|
||||||
substep_polyinj, substep_polyprod);
|
|
||||||
injected[0] += substep_injected[0];
|
|
||||||
injected[1] += substep_injected[1];
|
|
||||||
produced[0] += substep_produced[0];
|
|
||||||
produced[1] += substep_produced[1];
|
|
||||||
polyinj += substep_polyinj;
|
|
||||||
polyprod += substep_polyprod;
|
|
||||||
if (gravity_ != 0 && use_segregation_split_) {
|
|
||||||
tsolver_.solveGravity(columns_, stepsize,
|
|
||||||
state.saturation(), state.surfacevol(),
|
|
||||||
state.concentration(), state.maxconcentration());
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
transport_timer.stop();
|
} while (!well_control_passed);
|
||||||
double tt = transport_timer.secsSinceStart();
|
|
||||||
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
|
|
||||||
ttime += tt;
|
|
||||||
|
|
||||||
// Report volume balances.
|
// Update pore volumes if rock is compressible.
|
||||||
Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
|
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||||
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
|
initial_porevol = porevol;
|
||||||
polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
|
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
||||||
state, rock_comp_props_);
|
}
|
||||||
tot_injected[0] += injected[0];
|
|
||||||
tot_injected[1] += injected[1];
|
|
||||||
tot_produced[0] += produced[0];
|
|
||||||
tot_produced[1] += produced[1];
|
|
||||||
tot_polyinj += polyinj;
|
|
||||||
tot_polyprod += polyprod;
|
|
||||||
std::cout.precision(5);
|
|
||||||
const int width = 18;
|
|
||||||
std::cout << "\nMass balance: "
|
|
||||||
" water(surfvol) oil(surfvol) polymer(kg)\n";
|
|
||||||
std::cout << " In-place: "
|
|
||||||
<< std::setw(width) << inplace_surfvol[0]
|
|
||||||
<< std::setw(width) << inplace_surfvol[1]
|
|
||||||
<< std::setw(width) << polymass << std::endl;
|
|
||||||
std::cout << " Adsorbed: "
|
|
||||||
<< std::setw(width) << 0.0
|
|
||||||
<< std::setw(width) << 0.0
|
|
||||||
<< std::setw(width) << polymass_adsorbed << std::endl;
|
|
||||||
std::cout << " Injected: "
|
|
||||||
<< std::setw(width) << injected[0]
|
|
||||||
<< std::setw(width) << injected[1]
|
|
||||||
<< std::setw(width) << polyinj << std::endl;
|
|
||||||
std::cout << " Produced: "
|
|
||||||
<< std::setw(width) << produced[0]
|
|
||||||
<< std::setw(width) << produced[1]
|
|
||||||
<< std::setw(width) << polyprod << std::endl;
|
|
||||||
std::cout << " Total inj: "
|
|
||||||
<< std::setw(width) << tot_injected[0]
|
|
||||||
<< std::setw(width) << tot_injected[1]
|
|
||||||
<< std::setw(width) << tot_polyinj << std::endl;
|
|
||||||
std::cout << " Total prod: "
|
|
||||||
<< std::setw(width) << tot_produced[0]
|
|
||||||
<< std::setw(width) << tot_produced[1]
|
|
||||||
<< std::setw(width) << tot_polyprod << std::endl;
|
|
||||||
const double balance[3] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
|
|
||||||
init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1],
|
|
||||||
init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed };
|
|
||||||
std::cout << " Initial - inplace + inj - prod: "
|
|
||||||
<< std::setw(width) << balance[0]
|
|
||||||
<< std::setw(width) << balance[1]
|
|
||||||
<< std::setw(width) << balance[2]
|
|
||||||
<< std::endl;
|
|
||||||
std::cout << " Relative mass error: "
|
|
||||||
<< std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
|
|
||||||
<< std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
|
|
||||||
<< std::setw(width) << balance[2]/(init_polymass + tot_polyinj)
|
|
||||||
<< std::endl;
|
|
||||||
std::cout.precision(8);
|
|
||||||
|
|
||||||
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
|
// Process transport sources (to include bdy terms and well flows).
|
||||||
produced[0]/(produced[0] + produced[1]),
|
Opm::computeTransportSource(props_, wells_, well_state, transport_src);
|
||||||
tot_produced[0]/tot_porevol_init);
|
|
||||||
if (wells_) {
|
// Find inflow rate.
|
||||||
wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
|
const double current_time = timer.simulationTimeElapsed();
|
||||||
state.saturation(), timer.simulationTimeElapsed() + timer.currentStepLength(),
|
double stepsize = timer.currentStepLength();
|
||||||
well_state.bhp(), well_state.perfRates());
|
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
|
||||||
|
|
||||||
|
|
||||||
|
// Solve transport.
|
||||||
|
transport_timer.start();
|
||||||
|
if (num_transport_substeps_ != 1) {
|
||||||
|
stepsize /= double(num_transport_substeps_);
|
||||||
|
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
|
||||||
|
}
|
||||||
|
double injected[2] = { 0.0 };
|
||||||
|
double produced[2] = { 0.0 };
|
||||||
|
double polyinj = 0.0;
|
||||||
|
double polyprod = 0.0;
|
||||||
|
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
|
||||||
|
tsolver_.solve(&state.faceflux()[0], initial_pressure,
|
||||||
|
state.pressure(), &initial_porevol[0], &porevol[0],
|
||||||
|
&transport_src[0], &polymer_inflow_c[0], stepsize,
|
||||||
|
state.saturation(), state.surfacevol(),
|
||||||
|
state.concentration(), state.maxconcentration());
|
||||||
|
double substep_injected[2] = { 0.0 };
|
||||||
|
double substep_produced[2] = { 0.0 };
|
||||||
|
double substep_polyinj = 0.0;
|
||||||
|
double substep_polyprod = 0.0;
|
||||||
|
Opm::computeInjectedProduced(props_, poly_props_,
|
||||||
|
state,
|
||||||
|
transport_src, polymer_inflow_c, stepsize,
|
||||||
|
substep_injected, substep_produced,
|
||||||
|
substep_polyinj, substep_polyprod);
|
||||||
|
injected[0] += substep_injected[0];
|
||||||
|
injected[1] += substep_injected[1];
|
||||||
|
produced[0] += substep_produced[0];
|
||||||
|
produced[1] += substep_produced[1];
|
||||||
|
polyinj += substep_polyinj;
|
||||||
|
polyprod += substep_polyprod;
|
||||||
|
if (gravity_ != 0 && use_segregation_split_) {
|
||||||
|
tsolver_.solveGravity(columns_, stepsize,
|
||||||
|
state.saturation(), state.surfacevol(),
|
||||||
|
state.concentration(), state.maxconcentration());
|
||||||
}
|
}
|
||||||
// }
|
}
|
||||||
|
transport_timer.stop();
|
||||||
|
double tt = transport_timer.secsSinceStart();
|
||||||
|
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
|
||||||
|
ttime += tt;
|
||||||
|
|
||||||
|
// Report volume balances.
|
||||||
|
Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol);
|
||||||
|
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
|
||||||
|
polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_,
|
||||||
|
state, rock_comp_props_);
|
||||||
|
tot_injected[0] += injected[0];
|
||||||
|
tot_injected[1] += injected[1];
|
||||||
|
tot_produced[0] += produced[0];
|
||||||
|
tot_produced[1] += produced[1];
|
||||||
|
tot_polyinj += polyinj;
|
||||||
|
tot_polyprod += polyprod;
|
||||||
|
std::cout.precision(5);
|
||||||
|
const int width = 18;
|
||||||
|
std::cout << "\nMass balance: "
|
||||||
|
" water(surfvol) oil(surfvol) polymer(kg)\n";
|
||||||
|
std::cout << " In-place: "
|
||||||
|
<< std::setw(width) << inplace_surfvol[0]
|
||||||
|
<< std::setw(width) << inplace_surfvol[1]
|
||||||
|
<< std::setw(width) << polymass << std::endl;
|
||||||
|
std::cout << " Adsorbed: "
|
||||||
|
<< std::setw(width) << 0.0
|
||||||
|
<< std::setw(width) << 0.0
|
||||||
|
<< std::setw(width) << polymass_adsorbed << std::endl;
|
||||||
|
std::cout << " Injected: "
|
||||||
|
<< std::setw(width) << injected[0]
|
||||||
|
<< std::setw(width) << injected[1]
|
||||||
|
<< std::setw(width) << polyinj << std::endl;
|
||||||
|
std::cout << " Produced: "
|
||||||
|
<< std::setw(width) << produced[0]
|
||||||
|
<< std::setw(width) << produced[1]
|
||||||
|
<< std::setw(width) << polyprod << std::endl;
|
||||||
|
std::cout << " Total inj: "
|
||||||
|
<< std::setw(width) << tot_injected[0]
|
||||||
|
<< std::setw(width) << tot_injected[1]
|
||||||
|
<< std::setw(width) << tot_polyinj << std::endl;
|
||||||
|
std::cout << " Total prod: "
|
||||||
|
<< std::setw(width) << tot_produced[0]
|
||||||
|
<< std::setw(width) << tot_produced[1]
|
||||||
|
<< std::setw(width) << tot_polyprod << std::endl;
|
||||||
|
const double balance[3] = { init_surfvol[0] - inplace_surfvol[0] - tot_produced[0] + tot_injected[0],
|
||||||
|
init_surfvol[1] - inplace_surfvol[1] - tot_produced[1] + tot_injected[1],
|
||||||
|
init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed };
|
||||||
|
std::cout << " Initial - inplace + inj - prod: "
|
||||||
|
<< std::setw(width) << balance[0]
|
||||||
|
<< std::setw(width) << balance[1]
|
||||||
|
<< std::setw(width) << balance[2]
|
||||||
|
<< std::endl;
|
||||||
|
std::cout << " Relative mass error: "
|
||||||
|
<< std::setw(width) << balance[0]/(init_surfvol[0] + tot_injected[0])
|
||||||
|
<< std::setw(width) << balance[1]/(init_surfvol[1] + tot_injected[1])
|
||||||
|
<< std::setw(width) << balance[2]/(init_polymass + tot_polyinj)
|
||||||
|
<< std::endl;
|
||||||
|
std::cout.precision(8);
|
||||||
|
|
||||||
|
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
|
||||||
|
produced[0]/(produced[0] + produced[1]),
|
||||||
|
tot_produced[0]/tot_porevol_init);
|
||||||
|
if (wells_) {
|
||||||
|
wellreport.push(props_, *wells_, state.pressure(), state.surfacevol(),
|
||||||
|
state.saturation(), timer.simulationTimeElapsed() + timer.currentStepLength(),
|
||||||
|
well_state.bhp(), well_state.perfRates());
|
||||||
|
}
|
||||||
|
|
||||||
if (output_) {
|
if (output_) {
|
||||||
if (output_vtk_) {
|
if (output_vtk_) {
|
||||||
|
|||||||
@ -315,196 +315,194 @@ namespace Opm
|
|||||||
well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
|
well_resflows_phase.resize((wells_->number_of_phases)*(wells_->number_of_wells), 0.0);
|
||||||
wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
|
wellreport.push(props_, *wells_, state.saturation(), 0.0, well_state.bhp(), well_state.perfRates());
|
||||||
}
|
}
|
||||||
for (; !timer.done(); ++timer) {
|
// Report timestep and (optionally) write state to disk.
|
||||||
// Report timestep and (optionally) write state to disk.
|
timer.report(std::cout);
|
||||||
timer.report(std::cout);
|
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
|
||||||
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
|
if (output_vtk_) {
|
||||||
if (output_vtk_) {
|
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
}
|
||||||
|
if (output_binary_) {
|
||||||
|
outputStateBinary(grid_, state, timer, output_dir_);
|
||||||
|
}
|
||||||
|
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Solve pressure.
|
||||||
|
if (check_well_controls_) {
|
||||||
|
computeFractionalFlow(props_, poly_props_, allcells_,
|
||||||
|
state.saturation(), state.concentration(), state.maxconcentration(),
|
||||||
|
fractional_flows);
|
||||||
|
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
|
||||||
|
}
|
||||||
|
bool well_control_passed = !check_well_controls_;
|
||||||
|
int well_control_iteration = 0;
|
||||||
|
do {
|
||||||
|
// Run solver.
|
||||||
|
pressure_timer.start();
|
||||||
|
std::vector<double> initial_pressure = state.pressure();
|
||||||
|
psolver_.solve(timer.currentStepLength(), state, well_state);
|
||||||
|
|
||||||
|
// Renormalize pressure if rock is incompressible, and
|
||||||
|
// there are no pressure conditions (bcs or wells).
|
||||||
|
// It is deemed sufficient for now to renormalize
|
||||||
|
// using geometric volume instead of pore volume.
|
||||||
|
if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive())
|
||||||
|
&& allNeumannBCs(bcs_) && allRateWells(wells_)) {
|
||||||
|
// Compute average pressures of previous and last
|
||||||
|
// step, and total volume.
|
||||||
|
double av_prev_press = 0.0;
|
||||||
|
double av_press = 0.0;
|
||||||
|
double tot_vol = 0.0;
|
||||||
|
const int num_cells = grid_.number_of_cells;
|
||||||
|
for (int cell = 0; cell < num_cells; ++cell) {
|
||||||
|
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
|
||||||
|
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
|
||||||
|
tot_vol += grid_.cell_volumes[cell];
|
||||||
}
|
}
|
||||||
if (output_binary_) {
|
// Renormalization constant
|
||||||
outputStateBinary(grid_, state, timer, output_dir_);
|
const double ren_const = (av_prev_press - av_press)/tot_vol;
|
||||||
|
for (int cell = 0; cell < num_cells; ++cell) {
|
||||||
|
state.pressure()[cell] += ren_const;
|
||||||
|
}
|
||||||
|
const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
|
||||||
|
for (int well = 0; well < num_wells; ++well) {
|
||||||
|
well_state.bhp()[well] += ren_const;
|
||||||
}
|
}
|
||||||
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Solve pressure.
|
// Stop timer and report.
|
||||||
|
pressure_timer.stop();
|
||||||
|
double pt = pressure_timer.secsSinceStart();
|
||||||
|
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
|
||||||
|
ptime += pt;
|
||||||
|
|
||||||
|
// Optionally, check if well controls are satisfied.
|
||||||
if (check_well_controls_) {
|
if (check_well_controls_) {
|
||||||
computeFractionalFlow(props_, poly_props_, allcells_,
|
Opm::computePhaseFlowRatesPerWell(*wells_,
|
||||||
state.saturation(), state.concentration(), state.maxconcentration(),
|
well_state.perfRates(),
|
||||||
fractional_flows);
|
fractional_flows,
|
||||||
wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
|
well_resflows_phase);
|
||||||
}
|
std::cout << "Checking well conditions." << std::endl;
|
||||||
bool well_control_passed = !check_well_controls_;
|
// For testing we set surface := reservoir
|
||||||
int well_control_iteration = 0;
|
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
|
||||||
do {
|
++well_control_iteration;
|
||||||
// Run solver.
|
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
|
||||||
pressure_timer.start();
|
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
|
||||||
std::vector<double> initial_pressure = state.pressure();
|
|
||||||
psolver_.solve(timer.currentStepLength(), state, well_state);
|
|
||||||
|
|
||||||
// Renormalize pressure if rock is incompressible, and
|
|
||||||
// there are no pressure conditions (bcs or wells).
|
|
||||||
// It is deemed sufficient for now to renormalize
|
|
||||||
// using geometric volume instead of pore volume.
|
|
||||||
if ((rock_comp_props_ == NULL || !rock_comp_props_->isActive())
|
|
||||||
&& allNeumannBCs(bcs_) && allRateWells(wells_)) {
|
|
||||||
// Compute average pressures of previous and last
|
|
||||||
// step, and total volume.
|
|
||||||
double av_prev_press = 0.0;
|
|
||||||
double av_press = 0.0;
|
|
||||||
double tot_vol = 0.0;
|
|
||||||
const int num_cells = grid_.number_of_cells;
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
av_prev_press += initial_pressure[cell]*grid_.cell_volumes[cell];
|
|
||||||
av_press += state.pressure()[cell]*grid_.cell_volumes[cell];
|
|
||||||
tot_vol += grid_.cell_volumes[cell];
|
|
||||||
}
|
|
||||||
// Renormalization constant
|
|
||||||
const double ren_const = (av_prev_press - av_press)/tot_vol;
|
|
||||||
for (int cell = 0; cell < num_cells; ++cell) {
|
|
||||||
state.pressure()[cell] += ren_const;
|
|
||||||
}
|
|
||||||
const int num_wells = (wells_ == NULL) ? 0 : wells_->number_of_wells;
|
|
||||||
for (int well = 0; well < num_wells; ++well) {
|
|
||||||
well_state.bhp()[well] += ren_const;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
|
if (!well_control_passed) {
|
||||||
// Stop timer and report.
|
std::cout << "Well controls not passed, solving again." << std::endl;
|
||||||
pressure_timer.stop();
|
} else {
|
||||||
double pt = pressure_timer.secsSinceStart();
|
std::cout << "Well conditions met." << std::endl;
|
||||||
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;
|
|
||||||
ptime += pt;
|
|
||||||
|
|
||||||
// Optionally, check if well controls are satisfied.
|
|
||||||
if (check_well_controls_) {
|
|
||||||
Opm::computePhaseFlowRatesPerWell(*wells_,
|
|
||||||
well_state.perfRates(),
|
|
||||||
fractional_flows,
|
|
||||||
well_resflows_phase);
|
|
||||||
std::cout << "Checking well conditions." << std::endl;
|
|
||||||
// For testing we set surface := reservoir
|
|
||||||
well_control_passed = wells_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
|
|
||||||
++well_control_iteration;
|
|
||||||
if (!well_control_passed && well_control_iteration > max_well_control_iterations_) {
|
|
||||||
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
|
|
||||||
}
|
|
||||||
if (!well_control_passed) {
|
|
||||||
std::cout << "Well controls not passed, solving again." << std::endl;
|
|
||||||
} else {
|
|
||||||
std::cout << "Well conditions met." << std::endl;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
} while (!well_control_passed);
|
|
||||||
|
|
||||||
// Update pore volumes if rock is compressible.
|
|
||||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
|
||||||
initial_porevol = porevol;
|
|
||||||
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Process transport sources (to include bdy terms and well flows).
|
|
||||||
Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
|
|
||||||
wells_, well_state.perfRates(), transport_src);
|
|
||||||
|
|
||||||
// Find inflow rate.
|
|
||||||
const double current_time = timer.simulationTimeElapsed();
|
|
||||||
double stepsize = timer.currentStepLength();
|
|
||||||
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
|
|
||||||
|
|
||||||
// Solve transport.
|
|
||||||
transport_timer.start();
|
|
||||||
if (num_transport_substeps_ != 1) {
|
|
||||||
stepsize /= double(num_transport_substeps_);
|
|
||||||
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
|
|
||||||
}
|
|
||||||
double substep_injected[2] = { 0.0 };
|
|
||||||
double substep_produced[2] = { 0.0 };
|
|
||||||
double substep_polyinj = 0.0;
|
|
||||||
double substep_polyprod = 0.0;
|
|
||||||
injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
|
|
||||||
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
|
|
||||||
tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
|
|
||||||
state.saturation(), state.concentration(), state.maxconcentration());
|
|
||||||
Opm::computeInjectedProduced(props_, poly_props_,
|
|
||||||
state,
|
|
||||||
transport_src, polymer_inflow_c, stepsize,
|
|
||||||
substep_injected, substep_produced, substep_polyinj, substep_polyprod);
|
|
||||||
injected[0] += substep_injected[0];
|
|
||||||
injected[1] += substep_injected[1];
|
|
||||||
produced[0] += substep_produced[0];
|
|
||||||
produced[1] += substep_produced[1];
|
|
||||||
polyinj += substep_polyinj;
|
|
||||||
polyprod += substep_polyprod;
|
|
||||||
if (use_segregation_split_) {
|
|
||||||
tsolver_.solveGravity(columns_, &porevol[0], stepsize,
|
|
||||||
state.saturation(), state.concentration(), state.maxconcentration());
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
transport_timer.stop();
|
} while (!well_control_passed);
|
||||||
double tt = transport_timer.secsSinceStart();
|
|
||||||
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
|
|
||||||
ttime += tt;
|
|
||||||
|
|
||||||
// Report volume balances.
|
// Update pore volumes if rock is compressible.
|
||||||
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
|
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||||
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
|
initial_porevol = porevol;
|
||||||
polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
|
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
||||||
tot_injected[0] += injected[0];
|
}
|
||||||
tot_injected[1] += injected[1];
|
|
||||||
tot_produced[0] += produced[0];
|
|
||||||
tot_produced[1] += produced[1];
|
|
||||||
tot_polyinj += polyinj;
|
|
||||||
tot_polyprod += polyprod;
|
|
||||||
std::cout.precision(5);
|
|
||||||
const int width = 18;
|
|
||||||
std::cout << "\nVolume and polymer mass balance: "
|
|
||||||
" water(pv) oil(pv) polymer(kg)\n";
|
|
||||||
std::cout << " Saturated volumes: "
|
|
||||||
<< std::setw(width) << satvol[0]/tot_porevol_init
|
|
||||||
<< std::setw(width) << satvol[1]/tot_porevol_init
|
|
||||||
<< std::setw(width) << polymass << std::endl;
|
|
||||||
std::cout << " Adsorbed volumes: "
|
|
||||||
<< std::setw(width) << 0.0
|
|
||||||
<< std::setw(width) << 0.0
|
|
||||||
<< std::setw(width) << polymass_adsorbed << std::endl;
|
|
||||||
std::cout << " Injected volumes: "
|
|
||||||
<< std::setw(width) << injected[0]/tot_porevol_init
|
|
||||||
<< std::setw(width) << injected[1]/tot_porevol_init
|
|
||||||
<< std::setw(width) << polyinj << std::endl;
|
|
||||||
std::cout << " Produced volumes: "
|
|
||||||
<< std::setw(width) << produced[0]/tot_porevol_init
|
|
||||||
<< std::setw(width) << produced[1]/tot_porevol_init
|
|
||||||
<< std::setw(width) << polyprod << std::endl;
|
|
||||||
std::cout << " Total inj volumes: "
|
|
||||||
<< std::setw(width) << tot_injected[0]/tot_porevol_init
|
|
||||||
<< std::setw(width) << tot_injected[1]/tot_porevol_init
|
|
||||||
<< std::setw(width) << tot_polyinj << std::endl;
|
|
||||||
std::cout << " Total prod volumes: "
|
|
||||||
<< std::setw(width) << tot_produced[0]/tot_porevol_init
|
|
||||||
<< std::setw(width) << tot_produced[1]/tot_porevol_init
|
|
||||||
<< std::setw(width) << tot_polyprod << std::endl;
|
|
||||||
std::cout << " In-place + prod - inj: "
|
|
||||||
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
|
|
||||||
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init
|
|
||||||
<< std::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl;
|
|
||||||
std::cout << " Init - now - pr + inj: "
|
|
||||||
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
|
|
||||||
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
|
|
||||||
<< std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed)
|
|
||||||
<< std::endl;
|
|
||||||
std::cout.precision(8);
|
|
||||||
|
|
||||||
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
|
// Process transport sources (to include bdy terms and well flows).
|
||||||
produced[0]/(produced[0] + produced[1]),
|
Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
|
||||||
tot_produced[0]/tot_porevol_init);
|
wells_, well_state.perfRates(), transport_src);
|
||||||
if (wells_) {
|
|
||||||
wellreport.push(props_, *wells_, state.saturation(),
|
// Find inflow rate.
|
||||||
timer.simulationTimeElapsed() + timer.currentStepLength(),
|
const double current_time = timer.simulationTimeElapsed();
|
||||||
well_state.bhp(), well_state.perfRates());
|
double stepsize = timer.currentStepLength();
|
||||||
|
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
|
||||||
|
|
||||||
|
// Solve transport.
|
||||||
|
transport_timer.start();
|
||||||
|
if (num_transport_substeps_ != 1) {
|
||||||
|
stepsize /= double(num_transport_substeps_);
|
||||||
|
std::cout << "Making " << num_transport_substeps_ << " transport substeps." << std::endl;
|
||||||
|
}
|
||||||
|
double substep_injected[2] = { 0.0 };
|
||||||
|
double substep_produced[2] = { 0.0 };
|
||||||
|
double substep_polyinj = 0.0;
|
||||||
|
double substep_polyprod = 0.0;
|
||||||
|
injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0;
|
||||||
|
for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
|
||||||
|
tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize,
|
||||||
|
state.saturation(), state.concentration(), state.maxconcentration());
|
||||||
|
Opm::computeInjectedProduced(props_, poly_props_,
|
||||||
|
state,
|
||||||
|
transport_src, polymer_inflow_c, stepsize,
|
||||||
|
substep_injected, substep_produced, substep_polyinj, substep_polyprod);
|
||||||
|
injected[0] += substep_injected[0];
|
||||||
|
injected[1] += substep_injected[1];
|
||||||
|
produced[0] += substep_produced[0];
|
||||||
|
produced[1] += substep_produced[1];
|
||||||
|
polyinj += substep_polyinj;
|
||||||
|
polyprod += substep_polyprod;
|
||||||
|
if (use_segregation_split_) {
|
||||||
|
tsolver_.solveGravity(columns_, &porevol[0], stepsize,
|
||||||
|
state.saturation(), state.concentration(), state.maxconcentration());
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
transport_timer.stop();
|
||||||
|
double tt = transport_timer.secsSinceStart();
|
||||||
|
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;
|
||||||
|
ttime += tt;
|
||||||
|
|
||||||
|
// Report volume balances.
|
||||||
|
Opm::computeSaturatedVol(porevol, state.saturation(), satvol);
|
||||||
|
polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol());
|
||||||
|
polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration());
|
||||||
|
tot_injected[0] += injected[0];
|
||||||
|
tot_injected[1] += injected[1];
|
||||||
|
tot_produced[0] += produced[0];
|
||||||
|
tot_produced[1] += produced[1];
|
||||||
|
tot_polyinj += polyinj;
|
||||||
|
tot_polyprod += polyprod;
|
||||||
|
std::cout.precision(5);
|
||||||
|
const int width = 18;
|
||||||
|
std::cout << "\nVolume and polymer mass balance: "
|
||||||
|
" water(pv) oil(pv) polymer(kg)\n";
|
||||||
|
std::cout << " Saturated volumes: "
|
||||||
|
<< std::setw(width) << satvol[0]/tot_porevol_init
|
||||||
|
<< std::setw(width) << satvol[1]/tot_porevol_init
|
||||||
|
<< std::setw(width) << polymass << std::endl;
|
||||||
|
std::cout << " Adsorbed volumes: "
|
||||||
|
<< std::setw(width) << 0.0
|
||||||
|
<< std::setw(width) << 0.0
|
||||||
|
<< std::setw(width) << polymass_adsorbed << std::endl;
|
||||||
|
std::cout << " Injected volumes: "
|
||||||
|
<< std::setw(width) << injected[0]/tot_porevol_init
|
||||||
|
<< std::setw(width) << injected[1]/tot_porevol_init
|
||||||
|
<< std::setw(width) << polyinj << std::endl;
|
||||||
|
std::cout << " Produced volumes: "
|
||||||
|
<< std::setw(width) << produced[0]/tot_porevol_init
|
||||||
|
<< std::setw(width) << produced[1]/tot_porevol_init
|
||||||
|
<< std::setw(width) << polyprod << std::endl;
|
||||||
|
std::cout << " Total inj volumes: "
|
||||||
|
<< std::setw(width) << tot_injected[0]/tot_porevol_init
|
||||||
|
<< std::setw(width) << tot_injected[1]/tot_porevol_init
|
||||||
|
<< std::setw(width) << tot_polyinj << std::endl;
|
||||||
|
std::cout << " Total prod volumes: "
|
||||||
|
<< std::setw(width) << tot_produced[0]/tot_porevol_init
|
||||||
|
<< std::setw(width) << tot_produced[1]/tot_porevol_init
|
||||||
|
<< std::setw(width) << tot_polyprod << std::endl;
|
||||||
|
std::cout << " In-place + prod - inj: "
|
||||||
|
<< std::setw(width) << (satvol[0] + tot_produced[0] - tot_injected[0])/tot_porevol_init
|
||||||
|
<< std::setw(width) << (satvol[1] + tot_produced[1] - tot_injected[1])/tot_porevol_init
|
||||||
|
<< std::setw(width) << (polymass + tot_polyprod - tot_polyinj + polymass_adsorbed) << std::endl;
|
||||||
|
std::cout << " Init - now - pr + inj: "
|
||||||
|
<< std::setw(width) << (init_satvol[0] - satvol[0] - tot_produced[0] + tot_injected[0])/tot_porevol_init
|
||||||
|
<< std::setw(width) << (init_satvol[1] - satvol[1] - tot_produced[1] + tot_injected[1])/tot_porevol_init
|
||||||
|
<< std::setw(width) << (init_polymass - polymass - tot_polyprod + tot_polyinj - polymass_adsorbed)
|
||||||
|
<< std::endl;
|
||||||
|
std::cout.precision(8);
|
||||||
|
|
||||||
|
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
|
||||||
|
produced[0]/(produced[0] + produced[1]),
|
||||||
|
tot_produced[0]/tot_porevol_init);
|
||||||
|
if (wells_) {
|
||||||
|
wellreport.push(props_, *wells_, state.saturation(),
|
||||||
|
timer.simulationTimeElapsed() + timer.currentStepLength(),
|
||||||
|
well_state.bhp(), well_state.perfRates());
|
||||||
|
}
|
||||||
|
|
||||||
if (output_) {
|
if (output_) {
|
||||||
if (output_vtk_) {
|
if (output_vtk_) {
|
||||||
|
|||||||
@ -267,122 +267,105 @@ namespace Opm
|
|||||||
std::string filename = output_dir_ + "/step_timing.param";
|
std::string filename = output_dir_ + "/step_timing.param";
|
||||||
tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
|
tstep_os.open(filename.c_str(), std::fstream::out | std::fstream::app);
|
||||||
}
|
}
|
||||||
// while (!timer.done()) {
|
|
||||||
// Report timestep and (optionally) write state to disk.
|
// Report timestep and (optionally) write state to disk.
|
||||||
step_timer.start();
|
step_timer.start();
|
||||||
timer.report(std::cout);
|
timer.report(std::cout);
|
||||||
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
|
if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
|
||||||
if (output_vtk_) {
|
if (output_vtk_) {
|
||||||
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
}
|
|
||||||
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
|
||||||
// outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_);
|
|
||||||
|
|
||||||
}
|
}
|
||||||
|
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
|
}
|
||||||
|
|
||||||
SimulatorReport sreport;
|
SimulatorReport sreport;
|
||||||
|
|
||||||
// Solve pressure equation.
|
bool well_control_passed = !check_well_controls_;
|
||||||
// if (check_well_controls_) {
|
int well_control_iteration = 0;
|
||||||
// computeFractionalFlow(props_, allcells_,
|
do {
|
||||||
// state.pressure(), state.surfacevol(), state.saturation(),
|
// Run solver.
|
||||||
// fractional_flows);
|
const double current_time = timer.simulationTimeElapsed();
|
||||||
// wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase);
|
double stepsize = timer.currentStepLength();
|
||||||
// }
|
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
|
||||||
bool well_control_passed = !check_well_controls_;
|
solver_timer.start();
|
||||||
int well_control_iteration = 0;
|
std::vector<double> initial_pressure = state.pressure();
|
||||||
do {
|
solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src);
|
||||||
// Process transport sources (to include bdy terms and well flows).
|
// Stop timer and report.
|
||||||
// Opm::computeTransportSource(props_, wells_, well_state, transport_src);
|
solver_timer.stop();
|
||||||
// Run solver.
|
const double st = solver_timer.secsSinceStart();
|
||||||
const double current_time = timer.simulationTimeElapsed();
|
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
|
||||||
double stepsize = timer.currentStepLength();
|
|
||||||
polymer_inflow_.getInflowValues(current_time, current_time + stepsize, polymer_inflow_c);
|
|
||||||
solver_timer.start();
|
|
||||||
std::vector<double> initial_pressure = state.pressure();
|
|
||||||
solver_.step(timer.currentStepLength(), state, well_state, polymer_inflow_c, transport_src);
|
|
||||||
// Stop timer and report.
|
|
||||||
solver_timer.stop();
|
|
||||||
const double st = solver_timer.secsSinceStart();
|
|
||||||
std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl;
|
|
||||||
|
|
||||||
stime += st;
|
stime += st;
|
||||||
sreport.pressure_time = st;
|
sreport.pressure_time = st;
|
||||||
|
|
||||||
// Optionally, check if well controls are satisfied.
|
// Optionally, check if well controls are satisfied.
|
||||||
if (check_well_controls_) {
|
if (check_well_controls_) {
|
||||||
Opm::computePhaseFlowRatesPerWell(*wells_,
|
Opm::computePhaseFlowRatesPerWell(*wells_,
|
||||||
well_state.perfRates(),
|
well_state.perfRates(),
|
||||||
fractional_flows,
|
fractional_flows,
|
||||||
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_manager_.conditionsMet(well_state.bhp(), well_resflows_phase, well_resflows_phase);
|
well_control_passed = wells_manager_.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_) {
|
||||||
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
|
OPM_THROW(std::runtime_error, "Could not satisfy well conditions in " << max_well_control_iterations_ << " tries.");
|
||||||
}
|
|
||||||
if (!well_control_passed) {
|
|
||||||
std::cout << "Well controls not passed, solving again." << std::endl;
|
|
||||||
} else {
|
|
||||||
std::cout << "Well conditions met." << std::endl;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
} while (!well_control_passed);
|
if (!well_control_passed) {
|
||||||
// Update pore volumes if rock is compressible.
|
std::cout << "Well controls not passed, solving again." << std::endl;
|
||||||
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
} else {
|
||||||
initial_porevol = porevol;
|
std::cout << "Well conditions met." << std::endl;
|
||||||
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
|
||||||
}
|
|
||||||
|
|
||||||
double injected[2] = { 0.0 };
|
|
||||||
double produced[2] = { 0.0 };
|
|
||||||
double polyinj = 0;
|
|
||||||
double polyprod = 0;
|
|
||||||
|
|
||||||
Opm::computeInjectedProduced(props_, polymer_props_,
|
|
||||||
state,
|
|
||||||
transport_src, polymer_inflow_c, timer.currentStepLength(),
|
|
||||||
injected, produced,
|
|
||||||
polyinj, polyprod);
|
|
||||||
tot_injected[0] += injected[0];
|
|
||||||
tot_injected[1] += injected[1];
|
|
||||||
tot_produced[0] += produced[0];
|
|
||||||
tot_produced[1] += produced[1];
|
|
||||||
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
|
|
||||||
produced[0]/(produced[0] + produced[1]),
|
|
||||||
tot_produced[0]/tot_porevol_init);
|
|
||||||
std::cout.precision(5);
|
|
||||||
const int width = 18;
|
|
||||||
std::cout << "\nMass balance report.\n";
|
|
||||||
std::cout << " Injected reservoir volumes: "
|
|
||||||
<< std::setw(width) << injected[0]
|
|
||||||
<< std::setw(width) << injected[1] << std::endl;
|
|
||||||
std::cout << " Produced reservoir volumes: "
|
|
||||||
<< std::setw(width) << produced[0]
|
|
||||||
<< std::setw(width) << produced[1] << std::endl;
|
|
||||||
std::cout << " Total inj reservoir volumes: "
|
|
||||||
<< std::setw(width) << tot_injected[0]
|
|
||||||
<< std::setw(width) << tot_injected[1] << std::endl;
|
|
||||||
std::cout << " Total prod reservoir volumes: "
|
|
||||||
<< std::setw(width) << tot_produced[0]
|
|
||||||
<< std::setw(width) << tot_produced[1] << std::endl;
|
|
||||||
sreport.total_time = step_timer.secsSinceStart();
|
|
||||||
if (output_) {
|
|
||||||
sreport.reportParam(tstep_os);
|
|
||||||
|
|
||||||
if (output_vtk_) {
|
|
||||||
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
|
||||||
}
|
}
|
||||||
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
|
||||||
outputWaterCut(watercut, output_dir_);
|
|
||||||
tstep_os.close();
|
|
||||||
}
|
}
|
||||||
|
} while (!well_control_passed);
|
||||||
|
// Update pore volumes if rock is compressible.
|
||||||
|
if (rock_comp_props_ && rock_comp_props_->isActive()) {
|
||||||
|
initial_porevol = porevol;
|
||||||
|
computePorevolume(grid_, props_.porosity(), *rock_comp_props_, state.pressure(), porevol);
|
||||||
|
}
|
||||||
|
|
||||||
// advance to next timestep before reporting at this location
|
double injected[2] = { 0.0 };
|
||||||
// ++timer;
|
double produced[2] = { 0.0 };
|
||||||
|
double polyinj = 0;
|
||||||
|
double polyprod = 0;
|
||||||
|
|
||||||
// }
|
Opm::computeInjectedProduced(props_, polymer_props_,
|
||||||
|
state,
|
||||||
|
transport_src, polymer_inflow_c, timer.currentStepLength(),
|
||||||
|
injected, produced,
|
||||||
|
polyinj, polyprod);
|
||||||
|
tot_injected[0] += injected[0];
|
||||||
|
tot_injected[1] += injected[1];
|
||||||
|
tot_produced[0] += produced[0];
|
||||||
|
tot_produced[1] += produced[1];
|
||||||
|
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
|
||||||
|
produced[0]/(produced[0] + produced[1]),
|
||||||
|
tot_produced[0]/tot_porevol_init);
|
||||||
|
std::cout.precision(5);
|
||||||
|
const int width = 18;
|
||||||
|
std::cout << "\nMass balance report.\n";
|
||||||
|
std::cout << " Injected reservoir volumes: "
|
||||||
|
<< std::setw(width) << injected[0]
|
||||||
|
<< std::setw(width) << injected[1] << std::endl;
|
||||||
|
std::cout << " Produced reservoir volumes: "
|
||||||
|
<< std::setw(width) << produced[0]
|
||||||
|
<< std::setw(width) << produced[1] << std::endl;
|
||||||
|
std::cout << " Total inj reservoir volumes: "
|
||||||
|
<< std::setw(width) << tot_injected[0]
|
||||||
|
<< std::setw(width) << tot_injected[1] << std::endl;
|
||||||
|
std::cout << " Total prod reservoir volumes: "
|
||||||
|
<< std::setw(width) << tot_produced[0]
|
||||||
|
<< std::setw(width) << tot_produced[1] << std::endl;
|
||||||
|
sreport.total_time = step_timer.secsSinceStart();
|
||||||
|
if (output_) {
|
||||||
|
sreport.reportParam(tstep_os);
|
||||||
|
|
||||||
|
if (output_vtk_) {
|
||||||
|
outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
|
}
|
||||||
|
outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
|
||||||
|
outputWaterCut(watercut, output_dir_);
|
||||||
|
tstep_os.close();
|
||||||
|
}
|
||||||
|
|
||||||
total_timer.stop();
|
total_timer.stop();
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user