From 18246263e921516475a84f28749475dad0e2a210 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 12 Apr 2016 08:50:34 +0200 Subject: [PATCH 1/7] Move computation of well potentials from simulator class ot model class - the computation of well potentials in the model class calculates the well potentials using computeWellFlux() - in this way the well potential calculations also handle well where some perforations are closed by the simulator due to cross-flow. - the well potentials pr perforation and phase is stored in the well state. --- opm/autodiff/BlackoilModelBase.hpp | 22 +++++ opm/autodiff/BlackoilModelBase_impl.hpp | 95 +++++++++++++++++++ opm/autodiff/SimulatorBase_impl.hpp | 22 +---- .../WellStateFullyImplicitBlackoil.hpp | 8 ++ 4 files changed, 127 insertions(+), 20 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 41af2dc3c..4088e46f8 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -406,6 +406,28 @@ namespace Opm { WellState& well_state); void +<<<<<<< HEAD +======= + computeWellFlux(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + V& aliveWells, + std::vector& cq_s) const; + + void + updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const; + + void + computeWellPotentials(const SolutionState& state, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + WellState& well_state); + + + void +>>>>>>> Move computation of well potentials from simulator class ot model class 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..ef2922c6b 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -867,9 +867,104 @@ 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); + } + template + void + BlackoilModelBase:: + computeWellPotentials(const SolutionState& state, + const std::vector& mob_perfcells, + 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 largest potential 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]; + +// if (well_type == INJECTOR) { +// double dp = detail::computeHydrostaticCorrection( +// wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), +// well_perforation_densities_, gravity); +// const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; +// // pick the bhp that gives the largest potentials i.e. largest bhp for injectors +// if ( bhp > bhps[w]) { +// bhps[w] = bhp; +// } +// } +// else if (well_type == PRODUCER) { +// double dp = detail::computeHydrostaticCorrection( +// wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), +// well_perforation_densities_, gravity); + +// const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; +// // pick the bhp that gives the largest potentials i.e. smalest 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().computeWellFlux(state0, 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/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index c5d66d6d7..800f08a2e 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -396,30 +396,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 From cc48e8a17e96b3188d6386703c2921d8c55c2540 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 12 Apr 2016 09:19:55 +0200 Subject: [PATCH 2/7] Add support for THP in well potential calculations --- opm/autodiff/BlackoilModelBase_impl.hpp | 190 ++++++++++++------------ 1 file changed, 97 insertions(+), 93 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index ef2922c6b..df633f03c 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -872,99 +872,6 @@ namespace detail { } - template - void - BlackoilModelBase:: - computeWellPotentials(const SolutionState& state, - const std::vector& mob_perfcells, - 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 largest potential 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]; - -// if (well_type == INJECTOR) { -// double dp = detail::computeHydrostaticCorrection( -// wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), -// well_perforation_densities_, gravity); -// const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; -// // pick the bhp that gives the largest potentials i.e. largest bhp for injectors -// if ( bhp > bhps[w]) { -// bhps[w] = bhp; -// } -// } -// else if (well_type == PRODUCER) { -// double dp = detail::computeHydrostaticCorrection( -// wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), -// well_perforation_densities_, gravity); - -// const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; -// // pick the bhp that gives the largest potentials i.e. smalest 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().computeWellFlux(state0, 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); - } @@ -1789,6 +1696,103 @@ namespace detail { } + template + void + BlackoilModelBase:: + computeWellPotentials(const SolutionState& state, + const std::vector& mob_perfcells, + 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 largest potential 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 = detail::computeHydrostaticCorrection( + wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), + well_perforation_densities_, gravity); + const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + // pick the bhp that gives the largest potentials i.e. largest bhp for injectors + if ( bhp > bhps[w]) { + bhps[w] = bhp; + } + } + else if (well_type == PRODUCER) { + double dp = detail::computeHydrostaticCorrection( + wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), + well_perforation_densities_, gravity); + + const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + // pick the bhp that gives the largest potentials i.e. smalest 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().computeWellFlux(state0, 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); + } + + From b216302895bf9140928368c6b8be6769414973fb Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 12 Apr 2016 14:48:41 +0200 Subject: [PATCH 3/7] Adapt to usage of stdWells() --- opm/autodiff/BlackoilModelBase_impl.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index df633f03c..b80c1c13f 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -867,7 +867,6 @@ 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); } @@ -1747,7 +1746,7 @@ namespace detail { if (well_type == INJECTOR) { double dp = detail::computeHydrostaticCorrection( wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - well_perforation_densities_, gravity); + stdWells().wellPerforationDensities(), gravity); const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; // pick the bhp that gives the largest potentials i.e. largest bhp for injectors if ( bhp > bhps[w]) { @@ -1757,7 +1756,7 @@ namespace detail { else if (well_type == PRODUCER) { double dp = detail::computeHydrostaticCorrection( wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), - well_perforation_densities_, gravity); + stdWells().wellPerforationDensities(), gravity); const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; // pick the bhp that gives the largest potentials i.e. smalest bhp for producers From 4867e14f88f814b280126c2ff6a74ce56e4fc3af Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 19 Apr 2016 09:41:25 +0200 Subject: [PATCH 4/7] Fixes after rebase -- remove some leftovers from merge process -- start using stdWells() --- opm/autodiff/BlackoilModelBase.hpp | 15 --------------- opm/autodiff/BlackoilModelBase_impl.hpp | 6 +++--- 2 files changed, 3 insertions(+), 18 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 4088e46f8..9a59a7f64 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -406,20 +406,6 @@ namespace Opm { WellState& well_state); void -<<<<<<< HEAD -======= - computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s) const; - - void - updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - const SolutionState& state, - WellState& xw) const; - - void computeWellPotentials(const SolutionState& state, const std::vector& mob_perfcells, const std::vector& b_perfcells, @@ -427,7 +413,6 @@ namespace Opm { void ->>>>>>> Move computation of well potentials from simulator class ot model class addWellFluxEq(const std::vector& cq_s, const SolutionState& state); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index b80c1c13f..66dbb0e27 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1744,7 +1744,7 @@ namespace detail { const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); if (well_type == INJECTOR) { - double dp = detail::computeHydrostaticCorrection( + 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; @@ -1754,7 +1754,7 @@ namespace detail { } } else if (well_type == PRODUCER) { - double dp = detail::computeHydrostaticCorrection( + double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), stdWells().wellPerforationDensities(), gravity); @@ -1779,7 +1779,7 @@ namespace detail { // compute well potentials V aliveWells; std::vector well_potentials; - asImpl().computeWellFlux(state0, mob_perfcells, b_perfcells, aliveWells, 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 From 84972ac21e8f7d91939a78ec39d874618e6347d2 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 19 Apr 2016 15:58:05 +0200 Subject: [PATCH 5/7] Fix bhp/thp logic in well potential calculation Pick the most restrictive bhp (not the one that gives the largest potential flow) --- opm/autodiff/BlackoilModelBase_impl.hpp | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 66dbb0e27..ee76b7d7a 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1712,14 +1712,14 @@ namespace detail { 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 largest potential flow + //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) { + if (well_controls_iget_type(ctrl, ctrl_index) == THP) { double aqua = 0.0; double liquid = 0.0; double vapour = 0.0; @@ -1746,10 +1746,10 @@ namespace detail { if (well_type == INJECTOR) { double dp = wellhelpers::computeHydrostaticCorrection( wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); + stdWells().wellPerforationDensities(), gravity); const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - // pick the bhp that gives the largest potentials i.e. largest bhp for injectors - if ( bhp > bhps[w]) { + // apply the strictes of the bhp controlls i.e. smallest bhp for injectors + if ( bhp < bhps[w]) { bhps[w] = bhp; } } @@ -1759,8 +1759,8 @@ namespace detail { stdWells().wellPerforationDensities(), gravity); const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - // pick the bhp that gives the largest potentials i.e. smalest bhp for producers - if ( bhp < bhps[w]) { + // apply the strictes of the bhp controlls i.e. largest bhp for injectors + if ( bhp > bhps[w]) { bhps[w] = bhp; } } @@ -1769,6 +1769,7 @@ namespace detail { } } } + } // use bhp limit from control @@ -1789,6 +1790,7 @@ namespace detail { cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np); } well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np); + } From e9b097cbf4591898dec267d37993a310a66baaed Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 19 Apr 2016 16:03:33 +0200 Subject: [PATCH 6/7] 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. --- opm/autodiff/BlackoilModelBase_impl.hpp | 153 ++++++++++++----------- opm/autodiff/BlackoilModelParameters.cpp | 2 + opm/autodiff/BlackoilModelParameters.hpp | 4 + opm/autodiff/SimulatorBase_impl.hpp | 9 +- 4 files changed, 89 insertions(+), 79 deletions(-) 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 ); From 18434e54ad2686bb0fd01816b962c95e6d3b5d4e Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 20 Apr 2016 08:32:57 +0200 Subject: [PATCH 7/7] Query the compute_well_potential parameter only once -- and fix some comments --- opm/autodiff/BlackoilModelBase_impl.hpp | 4 ++-- opm/autodiff/SimulatorBase_impl.hpp | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 70164647f..dffd81769 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1750,7 +1750,7 @@ namespace detail { 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 + // apply the strictest of the bhp controlls i.e. smallest bhp for injectors if ( bhp < bhps[w]) { bhps[w] = bhp; } @@ -1761,7 +1761,7 @@ namespace detail { 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 + // apply the strictest of the bhp controlls i.e. largest bhp for producers if ( bhp > bhps[w]) { bhps[w] = bhp; } diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index c44c0025e..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,9 +222,9 @@ namespace Opm // Increment timer, remember well state. ++timer; prev_well_state = well_state; - // 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 ) ) { + // 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); }