diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 41af2dc3c..9a59a7f64 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -405,6 +405,13 @@ namespace Opm { SolutionState& state, WellState& well_state); + void + computeWellPotentials(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + WellState& well_state); + + void addWellFluxEq(const std::vector& cq_s, const SolutionState& state); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 7dc6dc2b8..dffd81769 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -867,8 +867,9 @@ namespace detail { asImpl().addWellFluxEq(cq_s, state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); asImpl().addWellControlEq(state, well_state, aliveWells); - } + asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state); + } @@ -1694,6 +1695,108 @@ namespace detail { } + template + void + BlackoilModelBase:: + computeWellPotentials(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + WellState& well_state) + { + //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 ] ]; + } + + 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); + + //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 strictest 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 strictest of the bhp controlls i.e. largest bhp for producers + 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); + } + + } + + 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 c5d66d6d7..915ad40e2 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -120,7 +120,7 @@ namespace Opm unsigned int totalNonlinearIterations = 0; unsigned int totalLinearIterations = 0; - + bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false ); std::vector well_potentials; // Main simulation loop. @@ -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); + // The well potentials are only computed if they are needed + // For now thay are only used to determine default guide rates for group controlled wells + if ( is_well_potentials_computed ) { + asImpl().computeWellPotentials(wells, state, well_state, well_potentials); + } } - // Write final simulation state. output_writer_.writeTimeStep( timer, state, prev_well_state ); @@ -396,30 +397,12 @@ namespace Opm { const int nw = wells->number_of_wells; const int np = wells->number_of_phases; - well_potentials.clear(); - well_potentials.resize(nw*np,0.0); + well_potentials.resize(nw*np,0.0); for (int w = 0; w < nw; ++w) { for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) { - const double well_cell_pressure = x.pressure()[wells->well_cells[perf]]; - const double drawdown_used = well_cell_pressure - xw.perfPress()[perf]; - const WellControls* ctrl = wells->ctrls[w]; - const int nwc = well_controls_get_num(ctrl); - //Loop over all controls until we find a BHP control - //that specifies what we need... - double bhp = 0.0; - for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { - if (well_controls_iget_type(ctrl, ctrl_index) == BHP) { - bhp = well_controls_iget_target(ctrl, ctrl_index); - } - // TODO: do something for thp; - } - // Calculate the pressure difference in the well perforation - const double dp = xw.perfPress()[perf] - xw.bhp()[w]; - const double drawdown_maximum = well_cell_pressure - (bhp + dp); - for (int phase = 0; phase < np; ++phase) { - well_potentials[w*np + phase] += (drawdown_maximum / drawdown_used * xw.perfPhaseRates()[perf*np + phase]); + well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase]; } } } diff --git a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp index 7912a106c..5cf888246 100644 --- a/opm/autodiff/WellStateFullyImplicitBlackoil.hpp +++ b/opm/autodiff/WellStateFullyImplicitBlackoil.hpp @@ -102,6 +102,9 @@ namespace Opm current_controls_[w] = well_controls_get_current(wells->ctrls[w]); } + well_potentials_.clear(); + well_potentials_.resize(nperf * np, 0.0); + // intialize wells that have been there before // order may change so the mapping is based on the well name if( ! prevState.wellMap().empty() ) @@ -184,9 +187,14 @@ namespace Opm std::vector& currentControls() { return current_controls_; } const std::vector& currentControls() const { return current_controls_; } + /// One rate per phase and well connection. + std::vector& wellPotentials() { return well_potentials_; } + const std::vector& wellPotentials() const { return well_potentials_; } + private: std::vector perfphaserates_; std::vector current_controls_; + std::vector well_potentials_; }; } // namespace Opm