Only compute well potential if needed

A boolen user parameter is added to controll the computation of well
potential.

This is a temporary fix to assure that no extra computation time is used
on well potential calculation if it is not needed. The long term fix
will require a more thorough revising of the well group implementation.
This commit is contained in:
Tor Harald Sandve
2016-04-19 16:03:33 +02:00
parent 84972ac21e
commit e9b097cbf4
4 changed files with 89 additions and 79 deletions

View File

@@ -1703,94 +1703,97 @@ namespace detail {
const std::vector<ADB>& b_perfcells,
WellState& well_state)
{
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
const Opm::PhaseUsage pu = fluid_.phaseUsage();
V bhps = V::Zero(nw);
for (int w = 0; w < nw; ++w) {
const WellControls* ctrl = wells().ctrls[w];
const int nwc = well_controls_get_num(ctrl);
//Loop over all controls until we find a BHP control
//or a THP control that specifies what we need.
//Pick the value that gives the most restrictive flow
for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
//only compute well potentials if they are needed
if (param_.compute_well_potentials_) {
const int nw = wells().number_of_wells;
const int np = wells().number_of_phases;
const Opm::PhaseUsage pu = fluid_.phaseUsage();
V bhps = V::Zero(nw);
for (int w = 0; w < nw; ++w) {
const WellControls* ctrl = wells().ctrls[w];
const int nwc = well_controls_get_num(ctrl);
//Loop over all controls until we find a BHP control
//or a THP control that specifies what we need.
//Pick the value that gives the most restrictive flow
for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
}
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
if (active_[ Water ]) {
aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
}
const int vfp = well_controls_iget_vfp(ctrl, ctrl_index);
const double& thp = well_controls_iget_target(ctrl, ctrl_index);
const double& alq = well_controls_iget_alq(ctrl, ctrl_index);
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
double aqua = 0.0;
double liquid = 0.0;
double vapour = 0.0;
//Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w];
if (active_[ Water ]) {
aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
}
if (active_[ Oil ]) {
liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
}
if (active_[ Gas ]) {
vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
}
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
const int vfp = well_controls_iget_vfp(ctrl, ctrl_index);
const double& thp = well_controls_iget_target(ctrl, ctrl_index);
const double& alq = well_controls_iget_alq(ctrl, ctrl_index);
if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
// apply the strictes of the bhp controlls i.e. smallest bhp for injectors
if ( bhp < bhps[w]) {
bhps[w] = bhp;
//Set *BHP* target by calculating bhp from THP
const WellType& well_type = wells().type[w];
const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
if (well_type == INJECTOR) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
// apply the strictes of the bhp controlls i.e. smallest bhp for injectors
if ( bhp < bhps[w]) {
bhps[w] = bhp;
}
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
// apply the strictes of the bhp controlls i.e. largest bhp for injectors
if ( bhp > bhps[w]) {
bhps[w] = bhp;
}
}
else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
else if (well_type == PRODUCER) {
double dp = wellhelpers::computeHydrostaticCorrection(
wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(),
stdWells().wellPerforationDensities(), gravity);
const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
// apply the strictes of the bhp controlls i.e. largest bhp for injectors
if ( bhp > bhps[w]) {
bhps[w] = bhp;
}
}
else {
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
}
}
}
// use bhp limit from control
SolutionState state0 = state;
asImpl().makeConstantState(state0);
state0.bhp = ADB::constant(bhps);
// compute well potentials
V aliveWells;
std::vector<ADB> well_potentials;
asImpl().stdWells().computeWellFlux(state0, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, well_potentials);
// store well potentials in the well state
// transform to a single vector instead of separate vectors pr phase
const int nperf = wells().well_connpos[nw];
V cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
}
well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
}
// use bhp limit from control
SolutionState state0 = state;
asImpl().makeConstantState(state0);
state0.bhp = ADB::constant(bhps);
// compute well potentials
V aliveWells;
std::vector<ADB> well_potentials;
asImpl().stdWells().computeWellFlux(state0, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, well_potentials);
// store well potentials in the well state
// transform to a single vector instead of separate vectors pr phase
const int nperf = wells().well_connpos[nw];
V cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
for (int phase = 1; phase < np; ++phase) {
cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
}
well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
}

View File

@@ -48,6 +48,7 @@ namespace Opm
tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_);
update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
compute_well_potentials_ = param.getDefault("compute_well_potentials", compute_well_potentials_);
}
@@ -65,6 +66,7 @@ namespace Opm
tolerance_wells_ = 1.0e-3;
solve_welleq_initially_ = true;
update_equations_scaling_ = false;
compute_well_potentials_ = false;
}

View File

@@ -49,6 +49,10 @@ namespace Opm
/// Update scaling factors for mass balance equations
bool update_equations_scaling_;
/// Compute well potentials, needed to calculate default guide rates for group
/// controlled wells
bool compute_well_potentials_;
/// Construct from user parameters or defaults.
explicit BlackoilModelParameters( const parameter::ParameterGroup& param );

View File

@@ -222,12 +222,13 @@ namespace Opm
// Increment timer, remember well state.
++timer;
prev_well_state = well_state;
// Compute Well potentials (only used to determine default guide rates for group controlled wells)
// TODO: add some logic to avoid unnecessary calulations of well potentials.
asImpl().computeWellPotentials(wells, state, well_state, well_potentials);
// Compute Well potentials if they are needed
// Only used to determine default guide rates for group controlled wells
if ( param_.getDefault("compute_well_potentials", false ) ) {
asImpl().computeWellPotentials(wells, state, well_state, well_potentials);
}
}
// Write final simulation state.
output_writer_.writeTimeStep( timer, state, prev_well_state );