diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index ee76b7d7a..70164647f 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1703,94 +1703,97 @@ namespace detail { const std::vector& 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 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 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); - } diff --git a/opm/autodiff/BlackoilModelParameters.cpp b/opm/autodiff/BlackoilModelParameters.cpp index 0cefa6139..5cdf1a2d5 100644 --- a/opm/autodiff/BlackoilModelParameters.cpp +++ b/opm/autodiff/BlackoilModelParameters.cpp @@ -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; } diff --git a/opm/autodiff/BlackoilModelParameters.hpp b/opm/autodiff/BlackoilModelParameters.hpp index 1d1e25525..f4a5dda74 100644 --- a/opm/autodiff/BlackoilModelParameters.hpp +++ b/opm/autodiff/BlackoilModelParameters.hpp @@ -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 ); diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index 800f08a2e..c44c0025e 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -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 );