From 393316600199ad1368b24d8b419d64a890ca8747 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 3 Oct 2023 13:28:23 +0200 Subject: [PATCH 1/5] Automatic choke --- .../wells/BlackoilWellModelGeneric.cpp | 4 +++- .../wells/BlackoilWellModel_impl.hpp | 8 +++++--- opm/simulators/wells/GroupState.cpp | 17 ++++++++++++++++ opm/simulators/wells/GroupState.hpp | 5 +++++ opm/simulators/wells/WellInterface_impl.hpp | 20 +++++++++++-------- 5 files changed, 42 insertions(+), 12 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 9df44e7a6..4eae4abc3 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1546,7 +1546,6 @@ updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor) } for (auto& well : well_container_generic_) { - // Producers only, since we so far only support the // "extended" network model (properties defined by // BRANPROP and NODEPROP) which only applies to producers. @@ -1563,6 +1562,9 @@ updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor) if (thp_is_limit) { ws.thp = well->getTHPConstraint(summaryState_); } + //PJPE: Set thp of wells belonging to a subsea manifold equal to the node_pressure + if (network.node(well->wellEcl().groupName()).as_choke()) + ws.thp = new_limit; } } } diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 89cf3ae38..79f9d95eb 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -46,6 +46,8 @@ #include #include +#include +#include #include #include @@ -1232,6 +1234,9 @@ namespace Opm { BlackoilWellModel:: updateWellControlsAndNetwork(const bool mandatory_network_balance, const double dt, DeferredLogger& local_deferredLogger) { + // PJPE: calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES) + computeWellGroupThp(local_deferredLogger); + // not necessarily that we always need to update once of the network solutions bool do_network_update = true; bool well_group_control_changed = false; @@ -2682,9 +2687,6 @@ namespace Opm { deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg); } } - // If we're using local well solves that include control switches, they also update - // operability, so reset before main iterations begin - well->resetWellOperability(); } updatePrimaryVariables(deferred_logger); diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index 008f772a0..262d33a72 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -105,6 +105,23 @@ GroupState::production_rates(const std::string& gname) const } //------------------------------------------------------------------------- +template +bool GroupState:: +GroupState::update_well_group_thp(const std::string& gname, const double& thp) +{ + this->group_thp[gname] = thp; +} + +template +double GroupState:: +GroupState::well_group_thp(const std::string& gname) const +{ + auto group_iter = this->group_thp.find(gname); + if (group_iter == this->group_thp.end()) + throw std::logic_error("No such group"); + + return group_iter->second; +} template void GroupState:: diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index 70c0bf61e..f3f01f70c 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -55,6 +55,7 @@ public: const std::vector& production_rates(const std::string& gname) const; void update_well_group_thp(const std::string& gname, const double& thp); + Scalar well_group_thp(const std::string& gname) const; bool has_production_reduction_rates(const std::string& gname) const; @@ -213,7 +214,11 @@ private: std::map inj_vrep_rate; std::map m_grat_sales_target; std::map m_gpmaint_target; +<<<<<<< HEAD std::map group_thp; +======= + std::map group_thp; +>>>>>>> 8e410ac73 (Automatic choke) std::map, Group::InjectionCMode> injection_controls; WellContainer gpmaint_state; diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 00983b534..0071749f7 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1473,17 +1473,21 @@ namespace Opm efficiencyFactor, deferred_logger); - // we don't want to scale with zero and get zero rates. - if (scale > 0) { - for (int p = 0; p 0) { + for (int p = 0; p Date: Thu, 2 Nov 2023 10:59:16 +0100 Subject: [PATCH 2/5] This is a combination of 3 commits. autochoke producers control mode set to THP simplyfying code as some logic is moved to the parser minor repair work --- .../wells/BlackoilWellModelGeneric.cpp | 3 --- .../wells/BlackoilWellModel_impl.hpp | 3 +++ opm/simulators/wells/WellInterface_impl.hpp | 19 ++++++++----------- 3 files changed, 11 insertions(+), 14 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 4eae4abc3..155118159 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1562,9 +1562,6 @@ updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor) if (thp_is_limit) { ws.thp = well->getTHPConstraint(summaryState_); } - //PJPE: Set thp of wells belonging to a subsea manifold equal to the node_pressure - if (network.node(well->wellEcl().groupName()).as_choke()) - ws.thp = new_limit; } } } diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 79f9d95eb..e78f01077 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -2686,6 +2686,9 @@ namespace Opm { const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the previous rates"; deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg); } + // If we're using local well solves that include control switches, they also update + // operability, so reset before main iterations begin + well->resetWellOperability(); } } updatePrimaryVariables(deferred_logger); diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 0071749f7..0b05c2a12 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1473,19 +1473,16 @@ namespace Opm efficiencyFactor, deferred_logger); - // we don't want to scale with zero and get zero rates. - if (scale > 0) { - for (int p = 0; p 0) { + for (int p = 0; p Date: Wed, 6 Dec 2023 17:31:51 +0100 Subject: [PATCH 3/5] allow individual well constraints before rebasing moved common thp calculation to updateWellControls Small fix clean up and improvements according reviewer comments clean up and improvements according reviewer comments, part 2 changed assessing safe THP range rebasing fixes removed unused argument rebasing --- .../wells/BlackoilWellModelGeneric.cpp | 1 + .../wells/BlackoilWellModel_impl.hpp | 38 ++++++++++++++++--- opm/simulators/wells/GroupState.cpp | 5 ++- opm/simulators/wells/WellBhpThpCalculator.hpp | 18 +++++++++ opm/simulators/wells/WellInterfaceGeneric.cpp | 7 ++++ opm/simulators/wells/WellInterface_impl.hpp | 4 +- 6 files changed, 64 insertions(+), 9 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 155118159..9df44e7a6 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1546,6 +1546,7 @@ updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor) } for (auto& well : well_container_generic_) { + // Producers only, since we so far only support the // "extended" network model (properties defined by // BRANPROP and NODEPROP) which only applies to producers. diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index e78f01077..747b7af07 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -55,6 +55,7 @@ #include #endif + #if HAVE_MPI #include #endif @@ -1234,9 +1235,6 @@ namespace Opm { BlackoilWellModel:: updateWellControlsAndNetwork(const bool mandatory_network_balance, const double dt, DeferredLogger& local_deferredLogger) { - // PJPE: calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES) - computeWellGroupThp(local_deferredLogger); - // not necessarily that we always need to update once of the network solutions bool do_network_update = true; bool well_group_control_changed = false; @@ -1331,6 +1329,7 @@ namespace Opm { const auto& balance = this->schedule()[reportStepIdx].network_balance(); const Scalar thp_tolerance = balance.thp_tolerance(); + if (!network.active()) { return; } @@ -1346,6 +1345,7 @@ namespace Opm { const auto ctrl = group.productionControls(summary_state); const auto cmode = ctrl.cmode; const auto pu = this->phase_usage_; + //TODO: Auto choke combined with RESV control is not supported std::vector resv_coeff(pu.num_phases, 1.0); Scalar gratTargetFromSales = 0.0; @@ -1376,6 +1376,31 @@ namespace Opm { return (group_rate - orig_target)/orig_target; }; + double min_thp, max_thp; + std::array range_initial; + //Find an initial bracket + if (!this->well_group_thp_calc_.has_value()){ + // Retrieve the terminal pressure of the associated root of the manifold group + std::string node_name = nodeName; + while (!network.node(node_name).terminal_pressure().has_value()) { + auto branch = network.uptree_branch(node_name).value(); + node_name = branch.uptree_node(); + } + + min_thp = network.node(node_name).terminal_pressure().value(); + std::optional approximate_solution0; + WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp, local_deferredLogger); + + // Narrow down the bracket + double low1, high1; + std::array range = {0.9*min_thp, 1.1*max_thp}; + std::optional appr_sol; + WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger); + min_thp = low1; + max_thp = high1; + range_initial = {min_thp, max_thp}; + } + const auto upbranch = network.uptree_branch(nodeName); const auto it = this->node_pressures_.find((*upbranch).uptree_node()); const Scalar nodal_pressure = it->second; @@ -1446,6 +1471,7 @@ namespace Opm { for (auto& well : this->well_container_) { std::string well_name = well->name(); if (group.hasWell(well_name)) { + well->setDynamicThpLimit(well_group_thp); } } @@ -2686,10 +2712,10 @@ namespace Opm { const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the previous rates"; deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg); } - // If we're using local well solves that include control switches, they also update - // operability, so reset before main iterations begin - well->resetWellOperability(); } + // If we're using local well solves that include control switches, they also update + // operability, so reset before main iterations begin + well->resetWellOperability(); } updatePrimaryVariables(deferred_logger); diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index 262d33a72..d73c7ed63 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -105,8 +105,9 @@ GroupState::production_rates(const std::string& gname) const } //------------------------------------------------------------------------- + template -bool GroupState:: +void GroupState:: GroupState::update_well_group_thp(const std::string& gname, const double& thp) { this->group_thp[gname] = thp; @@ -123,6 +124,8 @@ GroupState::well_group_thp(const std::string& gname) const return group_iter->second; } +//------------------------------------------------------------------------- + template void GroupState:: GroupState::update_well_group_thp(const std::string& gname, const double& thp) diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index 4ece9168a..a8ca0b735 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -130,6 +130,24 @@ public: static bool bruteForceBracketCommonTHP(const std::function& eq, Scalar& min_thp, Scalar& max_thp); + //! \brief Find limits using brute-force solver. + static bool bruteForceBracket(const std::function& eq, + const std::array& range, + double& low, double& high, + DeferredLogger& deferred_logger); + + //! \brief Find limits using brute-force solver. + static bool bruteForceBracketCommonTHP(const std::function& eq, + const std::array& range, + double& low, double& high, + std::optional& approximate_solution, + const double& limit, + DeferredLogger& deferred_logger); + + //! \brief Find limits using brute-force solver. + static bool bruteForceBracketCommonTHP(const std::function& eq, + double& min_thp, double& max_thp); + private: //! \brief Compute BHP from THP limit for an injector - implementation. template diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 8a9748768..07c2cb3b4 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -518,6 +518,13 @@ WellInterfaceGeneric::getDynamicThpLimit() const return dynamic_thp_limit_; } +template +void WellInterfaceGeneric:: +setDynamicThpLimit(const std::optional thp_limit) +{ + dynamic_thp_limit_ = thp_limit; +} + template void WellInterfaceGeneric:: updatePerforatedCell(std::vector& is_cell_perforated) diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 0b05c2a12..af95e9cd4 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1480,11 +1480,11 @@ namespace Opm } ws.trivial_target = false; } else { - ws.trivial_target = true; + ws.trivial_target = true; } break; - } + } case Well::ProducerCMode::CMODE_UNDEFINED: case Well::ProducerCMode::NONE: { From 9d735b8d6edb164a89c64d5eafc8a9d33d3ddc9c Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 13 May 2024 11:22:30 +0200 Subject: [PATCH 4/5] add group target calculation added temporary output mainly rebasing rebasing some further attempts fixed target calculations remove some case specific choices clean up some clean up generalised code for calculating target rate in groupControlledWells small rebase fix replaced by replaced by (2) --- opm/simulators/wells/BlackoilWellModel.hpp | 3 +- .../wells/BlackoilWellModelGeneric.cpp | 2 + .../wells/BlackoilWellModel_impl.hpp | 70 ++++++----- opm/simulators/wells/FractionCalculator.cpp | 6 +- opm/simulators/wells/FractionCalculator.hpp | 2 + opm/simulators/wells/GroupState.cpp | 27 ++-- opm/simulators/wells/GroupState.hpp | 6 +- opm/simulators/wells/WellBhpThpCalculator.cpp | 3 +- opm/simulators/wells/WellBhpThpCalculator.hpp | 18 --- opm/simulators/wells/WellGroupControls.cpp | 117 ++++++++++++++++++ opm/simulators/wells/WellGroupControls.hpp | 15 +++ opm/simulators/wells/WellGroupHelpers.cpp | 114 +++++++++++++++-- opm/simulators/wells/WellGroupHelpers.hpp | 8 ++ opm/simulators/wells/WellInterfaceGeneric.cpp | 2 +- 14 files changed, 299 insertions(+), 94 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index a047d28f2..99991c673 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -479,7 +479,8 @@ template class WellContributions; bool updateWellControlsAndNetwork(const bool mandatory_network_balance, const double dt, DeferredLogger& local_deferredLogger); - + + double computeWellGroupTarget(DeferredLogger& local_deferredLogger); void computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger); /// Update rank's notion of intersecting wells and their diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 9df44e7a6..7f0b368ee 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1313,6 +1313,7 @@ updateAndCommunicateGroupData(const int reportStepIdx, phase_usage_, guideRate_, well_state, + summaryState_, this->groupState(), groupTargetReduction); std::vector groupTargetReductionInj(numPhases(), 0.0); @@ -1323,6 +1324,7 @@ updateAndCommunicateGroupData(const int reportStepIdx, phase_usage_, guideRate_, well_state, + summaryState_, this->groupState(), groupTargetReductionInj); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 747b7af07..deb9a526e 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -42,12 +42,11 @@ #include #include #include +#include #include #include #include -#include -#include #include #include @@ -55,7 +54,6 @@ #include #endif - #if HAVE_MPI #include #endif @@ -1329,7 +1327,6 @@ namespace Opm { const auto& balance = this->schedule()[reportStepIdx].network_balance(); const Scalar thp_tolerance = balance.thp_tolerance(); - if (!network.active()) { return; } @@ -1342,20 +1339,51 @@ namespace Opm { if (has_choke) { const auto& summary_state = this->simulator_.vanguard().summaryState(); const Group& group = this->schedule().getGroup(nodeName, reportStepIdx); - const auto ctrl = group.productionControls(summary_state); - const auto cmode = ctrl.cmode; - const auto pu = this->phase_usage_; + const auto pu = this->phase_usage_; //TODO: Auto choke combined with RESV control is not supported std::vector resv_coeff(pu.num_phases, 1.0); Scalar gratTargetFromSales = 0.0; if (group_state.has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); + const auto ctrl = group.productionControls(summary_state); + auto cmode_tmp = ctrl.cmode; + Scalar target_tmp{0.0}; + bool fld_none = false; + if (cmode_tmp == Group::ProductionCMode::FLD || cmode_tmp == Group::ProductionCMode::NONE) { + fld_none = true; + // Target is set for an ancestor group. Target for autochoke group to be + // derived from via group guide rates + const Scalar efficiencyFactor = 1.0; + const Group& parentGroup = this->schedule().getGroup(group.parent(), reportStepIdx); + auto target = WellGroupControls::getAutoChokeGroupProductionTargetRate( + group.name(), + parentGroup, + well_state, + group_state, + this->schedule(), + summary_state, + resv_coeff, + efficiencyFactor, + reportStepIdx, + pu, + &this->guideRate_, + local_deferredLogger); + target_tmp = target.first; + cmode_tmp = target.second; + } + const auto cmode = cmode_tmp; WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff, gratTargetFromSales, nodeName, group_state, group.has_gpmaint_control(cmode)); - const Scalar orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); + if (!fld_none) + { + target_tmp = tcalc.groupTarget(ctrl, local_deferredLogger); + } + + const Scalar orig_target = target_tmp; + std::cout<< "Group: " << group.name() << " orig_target: " << orig_target*86400 << std::endl; auto mismatch = [&] (auto group_thp) { Scalar group_rate(0.0); @@ -1376,31 +1404,6 @@ namespace Opm { return (group_rate - orig_target)/orig_target; }; - double min_thp, max_thp; - std::array range_initial; - //Find an initial bracket - if (!this->well_group_thp_calc_.has_value()){ - // Retrieve the terminal pressure of the associated root of the manifold group - std::string node_name = nodeName; - while (!network.node(node_name).terminal_pressure().has_value()) { - auto branch = network.uptree_branch(node_name).value(); - node_name = branch.uptree_node(); - } - - min_thp = network.node(node_name).terminal_pressure().value(); - std::optional approximate_solution0; - WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp, local_deferredLogger); - - // Narrow down the bracket - double low1, high1; - std::array range = {0.9*min_thp, 1.1*max_thp}; - std::optional appr_sol; - WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger); - min_thp = low1; - max_thp = high1; - range_initial = {min_thp, max_thp}; - } - const auto upbranch = network.uptree_branch(nodeName); const auto it = this->node_pressures_.find((*upbranch).uptree_node()); const Scalar nodal_pressure = it->second; @@ -1471,7 +1474,6 @@ namespace Opm { for (auto& well : this->well_container_) { std::string well_name = well->name(); if (group.hasWell(well_name)) { - well->setDynamicThpLimit(well_group_thp); } } diff --git a/opm/simulators/wells/FractionCalculator.cpp b/opm/simulators/wells/FractionCalculator.cpp index 9e7ae0ff3..6af54ee58 100644 --- a/opm/simulators/wells/FractionCalculator.cpp +++ b/opm/simulators/wells/FractionCalculator.cpp @@ -37,6 +37,7 @@ FractionCalculator:: FractionCalculator(const Schedule& schedule, const WellState& well_state, const GroupState& group_state, + const SummaryState& summary_state, const int report_step, const GuideRate* guide_rate, const GuideRateModel::Target target, @@ -46,6 +47,7 @@ FractionCalculator(const Schedule& schedule, : schedule_(schedule) , well_state_(well_state) , group_state_(group_state) + , summary_state_(summary_state) , report_step_(report_step) , guide_rate_(guide_rate) , target_(target) @@ -123,6 +125,7 @@ guideRateSum(const Group& group, total_guide_rate += guideRate(child_group, always_included_child); } } + for (const std::string& child_well : group.wells()) { bool included = (child_well == always_included_child); if (is_producer_) { @@ -147,7 +150,7 @@ guideRate(const std::string& name, return WellGroupHelpers::getGuideRate(name, schedule_, well_state_, group_state_, report_step_, guide_rate_, target_, pu_); } else { - if (groupControlledWells(name, always_included_child) > 0) { + if ((groupControlledWells(name, always_included_child) > 0)) { if (is_producer_ && guide_rate_->has(name)) { return guide_rate_->get(name, target_, getGroupRateVector(name)); } else if (!is_producer_ && guide_rate_->has(name, injection_phase_)) { @@ -174,6 +177,7 @@ groupControlledWells(const std::string& group_name, return WellGroupHelpers::groupControlledWells(schedule_, well_state_, this->group_state_, + this->summary_state_, report_step_, group_name, always_included_child, diff --git a/opm/simulators/wells/FractionCalculator.hpp b/opm/simulators/wells/FractionCalculator.hpp index d421f5c3d..06962a393 100644 --- a/opm/simulators/wells/FractionCalculator.hpp +++ b/opm/simulators/wells/FractionCalculator.hpp @@ -41,6 +41,7 @@ public: FractionCalculator(const Schedule& schedule, const WellState& well_state, const GroupState& group_state, + const SummaryState& summary_state, const int report_step, const GuideRate* guide_rate, const GuideRateModel::Target target, @@ -65,6 +66,7 @@ private: const Schedule& schedule_; const WellState& well_state_; const GroupState& group_state_; + const SummaryState& summary_state_; int report_step_; const GuideRate* guide_rate_; GuideRateModel::Target target_; diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index d73c7ed63..f1a78796f 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -113,26 +113,6 @@ GroupState::update_well_group_thp(const std::string& gname, const double& thp) this->group_thp[gname] = thp; } -template -double GroupState:: -GroupState::well_group_thp(const std::string& gname) const -{ - auto group_iter = this->group_thp.find(gname); - if (group_iter == this->group_thp.end()) - throw std::logic_error("No such group"); - - return group_iter->second; -} - -//------------------------------------------------------------------------- - -template -void GroupState:: -GroupState::update_well_group_thp(const std::string& gname, const double& thp) -{ - this->group_thp[gname] = thp; -} - template Scalar GroupState:: GroupState::well_group_thp(const std::string& gname) const @@ -144,6 +124,13 @@ GroupState::well_group_thp(const std::string& gname) const return group_iter->second; } +template +bool GroupState:: +GroupState::is_autochoke_group(const std::string& gname) const +{ + return (this->group_thp.count(gname) > 0); +} + //------------------------------------------------------------------------- template diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index f3f01f70c..e98a77d9f 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -55,8 +55,8 @@ public: const std::vector& production_rates(const std::string& gname) const; void update_well_group_thp(const std::string& gname, const double& thp); - Scalar well_group_thp(const std::string& gname) const; + bool is_autochoke_group(const std::string& gname) const; bool has_production_reduction_rates(const std::string& gname) const; void update_production_reduction_rates(const std::string& gname, @@ -214,11 +214,7 @@ private: std::map inj_vrep_rate; std::map m_grat_sales_target; std::map m_gpmaint_target; -<<<<<<< HEAD std::map group_thp; -======= - std::map group_thp; ->>>>>>> 8e410ac73 (Automatic choke) std::map, Group::InjectionCMode> injection_controls; WellContainer gpmaint_state; diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 6ac4abe39..b60ece73c 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -1051,13 +1051,14 @@ bruteForceBracketCommonTHP(const std::function& eq, Scalar& min_thp, Scalar& max_thp) { bool bracket_found = false; - constexpr int sample_number = 1000; + constexpr int sample_number = 1000; constexpr Scalar interval = 1E5; Scalar eq_low = eq(min_thp); Scalar eq_high = 0.0; for (int i = 0; i < sample_number + 1; ++i) { max_thp = min_thp + interval * i; eq_high = eq(max_thp); + // std::cout << "max_thp: " << max_thp/1E5 << " eq_high: " << eq_high << " eq_low: " << eq_low << std::endl; if (eq_high * eq_low <= 0.) { bracket_found = true; min_thp = max_thp - interval; diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index a8ca0b735..4ece9168a 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -130,24 +130,6 @@ public: static bool bruteForceBracketCommonTHP(const std::function& eq, Scalar& min_thp, Scalar& max_thp); - //! \brief Find limits using brute-force solver. - static bool bruteForceBracket(const std::function& eq, - const std::array& range, - double& low, double& high, - DeferredLogger& deferred_logger); - - //! \brief Find limits using brute-force solver. - static bool bruteForceBracketCommonTHP(const std::function& eq, - const std::array& range, - double& low, double& high, - std::optional& approximate_solution, - const double& limit, - DeferredLogger& deferred_logger); - - //! \brief Find limits using brute-force solver. - static bool bruteForceBracketCommonTHP(const std::function& eq, - double& min_thp, double& max_thp); - private: //! \brief Compute BHP from THP limit for an injector - implementation. template diff --git a/opm/simulators/wells/WellGroupControls.cpp b/opm/simulators/wells/WellGroupControls.cpp index acb623083..a91e482bc 100644 --- a/opm/simulators/wells/WellGroupControls.cpp +++ b/opm/simulators/wells/WellGroupControls.cpp @@ -151,6 +151,7 @@ getGroupInjectionControl(const Group& group, WGHelpers::FractionCalculator fcalc(schedule, well_state, group_state, + summaryState, well_.currentStep(), well_.guideRate(), tcalc.guideTargetMode(), @@ -284,6 +285,7 @@ getGroupInjectionTargetRate(const Group& group, WGHelpers::FractionCalculator fcalc(schedule, well_state, group_state, + summaryState, well_.currentStep(), well_.guideRate(), tcalc.guideTargetMode(), @@ -401,6 +403,7 @@ getGroupProductionControl(const Group& group, WGHelpers::FractionCalculator fcalc(schedule, well_state, group_state, + summaryState, well_.currentStep(), well_.guideRate(), tcalc.guideTargetMode(), @@ -499,6 +502,7 @@ getGroupProductionTargetRate(const Group& group, WGHelpers::FractionCalculator fcalc(schedule, well_state, group_state, + summaryState, well_.currentStep(), well_.guideRate(), tcalc.guideTargetMode(), @@ -525,14 +529,17 @@ getGroupProductionTargetRate(const Group& group, // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. const std::size_t num_ancestors = chain.size() - 1; Scalar target = orig_target; + std::cout << "target: " << target*86400 << " time: " << well_.currentStep() << " orig" << std::endl; for (std::size_t ii = 0; ii < num_ancestors; ++ii) { if ((ii == 0) || well_.guideRate()->has(chain[ii])) { // Apply local reductions only at the control level // (top) and for levels where we have a specified // group guide rate. target -= localReduction(chain[ii]); + std::cout << "ii: " << ii << " group: " << chain[ii] << " LocRed:" << localReduction(chain[ii])*86400 << " target: " << target*86400 << std::endl; } target *= localFraction(chain[ii+1]); + std::cout << "ii: " << ii << " group: " << chain[ii+1] << " LocFrac: " << localFraction(chain[ii+1]) << " target: " << target*86400 << std::endl; } // Avoid negative target rates coming from too large local reductions. const Scalar target_rate = std::max(Scalar(0.0), target / efficiencyFactor); @@ -549,6 +556,116 @@ getGroupProductionTargetRate(const Group& group, return scale; } +template +std::pair WellGroupControls:: +getAutoChokeGroupProductionTargetRate(const std::string& name, + const Group& group, + const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const std::vector& resv_coeff, + Scalar efficiencyFactor, + const int reportStepIdx, + const PhaseUsage& pu, + const GuideRate* guideRate, + DeferredLogger& deferred_logger) +{ + const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name()); + if (currentGroupControl == Group::ProductionCMode::FLD || + currentGroupControl == Group::ProductionCMode::NONE) { + if (!group.productionGroupControlAvailable()) { + return std::make_pair(1.0, currentGroupControl); + } else { + // Produce share of parents control + const auto& parent = schedule.getGroup(group.parent(), reportStepIdx); + efficiencyFactor *= group.getGroupEfficiencyFactor(); + return getAutoChokeGroupProductionTargetRate(name, parent, well_state, group_state, + schedule, summaryState, + resv_coeff, efficiencyFactor, reportStepIdx, pu, + guideRate, deferred_logger); + } + } + + if (!group.isProductionGroup()) { + return std::make_pair(1.0, currentGroupControl); + } + + // If we are here, we are at the topmost group to be visited in the recursion. + // This is the group containing the control we will check against. + + // Make conversion factors for RESV <-> surface rates. + // std::vector resv_coeff(well_.phaseUsage().num_phases, 1.0); + // rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS. + + // gconsale may adjust the grat target. + // the adjusted rates is send to the targetCalculator + Scalar gratTargetFromSales = 0.0; + if (group_state.has_grat_sales_target(group.name())) + gratTargetFromSales = group_state.grat_sales_target(group.name()); + + WGHelpers::TargetCalculator tcalc(currentGroupControl, + pu, + resv_coeff, + gratTargetFromSales, + group.name(), + group_state, + group.has_gpmaint_control(currentGroupControl)); + + WGHelpers::FractionCalculator fcalc(schedule, + well_state, + group_state, + summaryState, + reportStepIdx, + guideRate, + tcalc.guideTargetMode(), + pu, + true, + Phase::OIL); + + auto localFraction = [&](const std::string& child) { + return fcalc.localFraction(child, child); //Note child needs to be passed to always include since the global isGrup map is not updated yet. + }; + + auto localReduction = [&](const std::string& group_name) { + const std::vector& groupTargetReductions = group_state.production_reduction_rates(group_name); + return tcalc.calcModeRateFromRates(groupTargetReductions); + }; + + std::optional ctrl; + if (!group.has_gpmaint_control(currentGroupControl)) + ctrl = group.productionControls(summaryState); + + // Scalar fr_true = fcalc.fraction("B1", "M5S", true); + // fr_true = fcalc.fraction("B1", "M5S", true); + // Scalar fr_false = fcalc.fraction("B1", "M5S", false); + // fr_false = fcalc.fraction("B1", "M5S", false); + // std::cout << "fr_true: " << fr_true << " fr_false: " << fr_false << std::endl; + + const double orig_target = tcalc.groupTarget(ctrl, deferred_logger); + const auto chain = WellGroupHelpers::groupChainTopBot(name, group.name(), + schedule, reportStepIdx); + // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. + const std::size_t num_ancestors = chain.size() - 1; + double target = orig_target; + std::cout << "target: " << target*86400 << " time: " << reportStepIdx << " modified" << std::endl; + for (std::size_t ii = 0; ii < num_ancestors; ++ii) { + if ((ii == 0) || guideRate->has(chain[ii])) { + // Apply local reductions only at the control level + // (top) and for levels where we have a specified + // group guide rate. + target -= localReduction(chain[ii]); + std::cout << "ii: " << ii << " group: " << chain[ii] << " LocRed:" << localReduction(chain[ii])*86400 << " target: " << target*86400 << std::endl; + } + target *= localFraction(chain[ii+1]); + std::cout << "ii: " << ii << " group: " << chain[ii+1] << " LocFrac: " << localFraction(chain[ii+1]) << " target: " << target*86400 << std::endl; + } + // Avoid negative target rates coming from too large local reductions. + const double target_rate = std::max(0.0, target / efficiencyFactor); + + return std::make_pair(target_rate, currentGroupControl); +} + #define INSTANTIATE(T,...) \ template void WellGroupControls:: \ getGroupInjectionControl(const Group&, \ diff --git a/opm/simulators/wells/WellGroupControls.hpp b/opm/simulators/wells/WellGroupControls.hpp index 77ca5e4b8..93ec7030d 100644 --- a/opm/simulators/wells/WellGroupControls.hpp +++ b/opm/simulators/wells/WellGroupControls.hpp @@ -24,6 +24,7 @@ #ifndef OPM_WELL_GROUP_CONTROLS_HEADER_INCLUDED #define OPM_WELL_GROUP_CONTROLS_HEADER_INCLUDED +#include #include #include #include @@ -37,6 +38,7 @@ class Group; template class GroupState; enum class InjectorType; using RegionId = int; +struct PhaseUsage; class Schedule; class SummaryState; template class WellInterfaceGeneric; @@ -99,6 +101,19 @@ public: Scalar efficiencyFactor, DeferredLogger& deferred_logger) const; + static std::pair getAutoChokeGroupProductionTargetRate(const std::string& name, + const Group& parent, + const WellState& well_state, + const GroupState& group_state, + const Schedule& schedule, + const SummaryState& summaryState, + const std::vector& resv_coeff, + Scalar efficiencyFactor, + const int reportStepIdx, + const PhaseUsage& pu, + const GuideRate* guideRate, + DeferredLogger& deferred_logger); + private: const WellInterfaceGeneric& well_; //!< Reference to well interface }; diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 89363162b..bb6a133ff 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -35,6 +36,7 @@ #include #include +#include #include #include #include @@ -345,6 +347,7 @@ updateGroupTargetReduction(const Group& group, const PhaseUsage& pu, const GuideRate& guide_rate, const WellState& wellState, + const SummaryState& summaryState, GroupState& group_state, std::vector& groupTargetReduction) { @@ -359,6 +362,7 @@ updateGroupTargetReduction(const Group& group, pu, guide_rate, wellState, + summaryState, group_state, subGroupTargetReduction); @@ -400,7 +404,7 @@ updateGroupTargetReduction(const Group& group, const bool individual_control = (currentGroupControl != Group::InjectionCMode::FLD && currentGroupControl != Group::InjectionCMode::NONE); const int num_group_controlled_wells - = groupControlledWells(schedule, wellState, group_state, reportStepIdx, subGroupName, "", !isInjector, phase); + = groupControlledWells(schedule, wellState, group_state, summaryState, reportStepIdx, subGroupName, "", !isInjector, phase); if (individual_control || num_group_controlled_wells == 0) { groupTargetReduction[phase_pos] += subGroupEfficiency * sumWellSurfaceRates(subGroup, schedule, wellState, reportStepIdx, phase_pos, isInjector); @@ -416,7 +420,8 @@ updateGroupTargetReduction(const Group& group, const bool individual_control = (currentGroupControl != Group::ProductionCMode::FLD && currentGroupControl != Group::ProductionCMode::NONE); const int num_group_controlled_wells - = groupControlledWells(schedule, wellState, group_state, reportStepIdx, subGroupName, "", !isInjector, /*injectionPhaseNotUsed*/Phase::OIL); + = groupControlledWells(schedule, wellState, group_state, summaryState, reportStepIdx, subGroupName, "", !isInjector, /*injectionPhaseNotUsed*/Phase::OIL); + std::cout<< "Group: " << subGroupName << " NumGroupControlledWells: " << num_group_controlled_wells << std::endl; if (individual_control || num_group_controlled_wells == 0) { for (int phase = 0; phase < np; phase++) { groupTargetReduction[phase] @@ -431,6 +436,7 @@ updateGroupTargetReduction(const Group& group, } } } + std::cout << "A: Group: " << group.name() << ", currentGroupControl: " << Group::ProductionCMode2String(currentGroupControl) << ", GroupTargetReduction[1]: " << groupTargetReduction[1]*86400 << std::endl; } } @@ -464,10 +470,16 @@ updateGroupTargetReduction(const Group& group, groupTargetReduction[phase] += ws.surface_rates[phase] * efficiency; } } else { - if (ws.production_cmode != Well::ProducerCMode::GRUP) - for (int phase = 0; phase < np; phase++) { - groupTargetReduction[phase] -= ws.surface_rates[phase] * efficiency; + if ((ws.production_cmode != Well::ProducerCMode::GRUP)){ + if (!group.as_choke()) { + for (int phase = 0; phase < np; phase++) { + groupTargetReduction[phase] -= ws.surface_rates[phase] * efficiency; + } } + std::cout<< "B: " << " Group: " << group.name() << ", Well: " + << wellName << ", cmode: " << ws.production_cmode << ", GroupTargetReduction[1]: " + << groupTargetReduction[1]*86400 << std::endl; + } } } if (isInjector) @@ -1086,6 +1098,7 @@ int WellGroupHelpers:: groupControlledWells(const Schedule& schedule, const WellState& well_state, const GroupState& group_state, + const SummaryState& summary_state, const int report_step, const std::string& group_name, const std::string& always_included_child, @@ -1107,16 +1120,65 @@ groupControlledWells(const Schedule& schedule, if (included) { num_wells - += groupControlledWells(schedule, well_state, group_state, report_step, child_group, always_included_child, is_production_group, injection_phase); + += groupControlledWells(schedule, well_state, group_state, summary_state, report_step, child_group, always_included_child, is_production_group, injection_phase); } } for (const std::string& child_well : group.wells()) { bool included = (child_well == always_included_child); if (is_production_group) { - included = included || well_state.isProductionGrup(child_well); + included = included || well_state.isProductionGrup(child_well) || group.as_choke(); } else { included = included || well_state.isInjectionGrup(child_well); } + const auto ctrl1 = group_state.production_control(group.name()); + if (group.as_choke() && ((ctrl1 == Group::ProductionCMode::FLD) || (ctrl1 == Group::ProductionCMode::NONE))){ + // The auto choke group has not own group control but inherits control from an ancestor group. + // Number of wells should be calculated as zero when wells of auto choke group do not deliver target. + // This behaviour is then similar to no-autochoke group with wells not on GRUP control. + // The rates of these wells are summed up. The parent group target is reduced with this rate. + // This reduced target becomes the target of the other child group of this parent. + const PhaseUsage& pu = well_state.phaseUsage(); + std::vector rates(pu.num_phases, 0.0); + for (int phase_pos = 0; phase_pos < pu.num_phases; ++phase_pos) { + rates[phase_pos] = WellGroupHelpers::sumWellSurfaceRates(group, + schedule, + well_state, + report_step, + phase_pos, + false); + } + + // Get the ancestor of the auto choke group that has group control (cmode != FLD, NONE) + const auto& control_group_name = control_group(group, group_state, report_step, schedule); + const auto& control_group = schedule.getGroup(control_group_name, report_step); + const auto& ctrl = control_group.productionControls(summary_state); + const auto& control_group_guide_rate = ctrl.guide_rate; + const auto& control_group_cmode = ctrl.cmode; + + const auto& group_guide_rate = group.productionControls(summary_state).guide_rate; + + Scalar gratTargetFromSales = 0.0; + if (group_state.has_grat_sales_target(control_group_name)) + gratTargetFromSales = group_state.grat_sales_target(control_group_name); + + std::vector resv_coeff(pu.num_phases, 1.0); + WGHelpers::TargetCalculator tcalc(control_group_cmode, + pu, + resv_coeff, + gratTargetFromSales, + group.name(), + group_state, + group.has_gpmaint_control(control_group_cmode)); + auto deferred_logger = Opm::DeferredLogger(); + const auto& control_group_target = tcalc.groupTarget(ctrl, deferred_logger); + // Target rate for the auto choke group + const Scalar target_rate = control_group_target * group_guide_rate / control_group_guide_rate; + const Scalar current_rate = tcalc.calcModeRateFromRates(rates); + + if (current_rate < 0.99 * target_rate) + included = false; + } + if (included) { ++num_wells; } @@ -1155,6 +1217,27 @@ groupChainTopBot(const std::string& bottom, return chain; } +template +std::string +WellGroupHelpers:: +control_group(const Group& group, + const GroupState& group_state, + const int reportStepIdx, + const Schedule& schedule) +{ + const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name()); + + if (currentGroupControl == Group::ProductionCMode::FLD || currentGroupControl == Group::ProductionCMode::NONE) { + const auto& parent = schedule.getGroup(group.parent(), reportStepIdx); + return control_group(parent, + group_state, + reportStepIdx, + schedule); + } + + return group.name(); +} + template std::pair WellGroupHelpers:: @@ -1234,6 +1317,7 @@ checkGroupConstraintsProd(const std::string& name, WGHelpers::FractionCalculator fcalc(schedule, wellState, group_state, + summaryState, reportStepIdx, guideRate, tcalc.guideTargetMode(), @@ -1276,6 +1360,7 @@ checkGroupConstraintsProd(const std::string& name, const int num_gr_ctrl = groupControlledWells(schedule, wellState, group_state, + summaryState, reportStepIdx, chain[ii], "", @@ -1296,6 +1381,7 @@ checkGroupConstraintsProd(const std::string& name, // we add a factor here to avoid switching due to numerical instability const Scalar factor = 1.01; if (currentRateFraction > (guiderateFraction * factor)) { + std::cout << guided_group << " localCurrentRate(guided_group): " << localCurrentRate(guided_group) << chain[ii-1] << " localCurrentRate(chain[ii-1]) " << localCurrentRate(chain[ii-1]) << std::endl; return std::make_pair(true, guiderateFraction / currentRateFraction); } } @@ -1409,6 +1495,7 @@ checkGroupConstraintsInj(const std::string& name, WGHelpers::FractionCalculator fcalc(schedule, wellState, group_state, + summaryState, reportStepIdx, guideRate, tcalc.guideTargetMode(), @@ -1451,12 +1538,13 @@ checkGroupConstraintsInj(const std::string& name, for (std::size_t ii = 1; ii < num_ancestors; ++ii) { const int num_gr_ctrl = groupControlledWells(schedule, wellState, - group_state, - reportStepIdx, - chain[ii], - "", - /*is_producer*/ false, - injectionPhase); + group_state, + summaryState, + reportStepIdx, + chain[ii], + "", + /*is_producer*/ false, + injectionPhase); if (guideRate->has(chain[ii], injectionPhase) && num_gr_ctrl > 0) { local_reduction_level = ii; } diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 4fc906a08..9687954b0 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -105,6 +105,7 @@ public: const PhaseUsage& pu, const GuideRate& guide_rate, const WellState& wellState, + const SummaryState& summaryState, GroupState& group_state, std::vector& groupTargetReduction); @@ -247,6 +248,7 @@ public: static int groupControlledWells(const Schedule& schedule, const WellState& well_state, const GroupState& group_state, + const SummaryState& summary_state, const int report_step, const std::string& group_name, const std::string& always_included_child, @@ -276,6 +278,12 @@ public: const Schedule& schedule, const int report_step); + static std::string + control_group(const Group& group, + const GroupState& group_state, + const int reportStepIdx, + const Schedule& schedule); + static std::pair checkGroupConstraintsProd(const std::string& name, const std::string& parent, diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 07c2cb3b4..26ffb2f36 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -519,7 +519,7 @@ WellInterfaceGeneric::getDynamicThpLimit() const } template -void WellInterfaceGeneric:: +void WellInterfaceGeneric:: setDynamicThpLimit(const std::optional thp_limit) { dynamic_thp_limit_ = thp_limit; From 8158ae8abe9da66ac1971494b169f5efcebd2324 Mon Sep 17 00:00:00 2001 From: plgbrts Date: Thu, 19 Dec 2024 11:55:03 +0000 Subject: [PATCH 5/5] allow defaulted group guide rates clean up --- opm/simulators/wells/BlackoilWellModel.hpp | 2 +- .../wells/BlackoilWellModel_impl.hpp | 24 ++++--- opm/simulators/wells/FractionCalculator.cpp | 3 +- opm/simulators/wells/WellBhpThpCalculator.cpp | 1 - opm/simulators/wells/WellGroupControls.cpp | 12 ---- opm/simulators/wells/WellGroupHelpers.cpp | 69 +++++++++++-------- opm/simulators/wells/WellGroupHelpers.hpp | 1 + opm/simulators/wells/WellInterface_impl.hpp | 1 - 8 files changed, 61 insertions(+), 52 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 99991c673..746626aaf 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -481,7 +481,7 @@ template class WellContributions; DeferredLogger& local_deferredLogger); double computeWellGroupTarget(DeferredLogger& local_deferredLogger); - void computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger); + bool computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger); /// Update rank's notion of intersecting wells and their /// associate solution variables. diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index deb9a526e..f113197e2 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1318,7 +1318,7 @@ namespace Opm { // The wells of such group should have a common THP and total phase rate(s) obeying (if possible) // the well group constraint set by GCONPROD template - void + bool BlackoilWellModel:: computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger) { @@ -1328,12 +1328,13 @@ namespace Opm { const Scalar thp_tolerance = balance.thp_tolerance(); if (!network.active()) { - return; + return false; } auto& well_state = this->wellState(); auto& group_state = this->groupState(); + bool well_group_thp_updated = false; for (const std::string& nodeName : network.node_names()) { const bool has_choke = network.node(nodeName).as_choke(); if (has_choke) { @@ -1354,7 +1355,7 @@ namespace Opm { if (cmode_tmp == Group::ProductionCMode::FLD || cmode_tmp == Group::ProductionCMode::NONE) { fld_none = true; // Target is set for an ancestor group. Target for autochoke group to be - // derived from via group guide rates + // derived via group guide rates const Scalar efficiencyFactor = 1.0; const Group& parentGroup = this->schedule().getGroup(group.parent(), reportStepIdx); auto target = WellGroupControls::getAutoChokeGroupProductionTargetRate( @@ -1379,11 +1380,11 @@ namespace Opm { group.has_gpmaint_control(cmode)); if (!fld_none) { + // Target is set for the autochoke group itself target_tmp = tcalc.groupTarget(ctrl, local_deferredLogger); } const Scalar orig_target = target_tmp; - std::cout<< "Group: " << group.name() << " orig_target: " << orig_target*86400 << std::endl; auto mismatch = [&] (auto group_thp) { Scalar group_rate(0.0); @@ -1479,9 +1480,14 @@ namespace Opm { } // Use the group THP in computeNetworkPressures(). - group_state.update_well_group_thp(nodeName, well_group_thp); + const auto& current_well_group_thp = group_state.is_autochoke_group(nodeName) ? group_state.well_group_thp(nodeName) : 1e30; + if (std::abs(current_well_group_thp - well_group_thp) > balance.pressure_tolerance()) { + well_group_thp_updated = true; + group_state.update_well_group_thp(nodeName, well_group_thp); + } } } + return well_group_thp_updated; } template @@ -2237,19 +2243,21 @@ namespace Opm { if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) { const double dt = this->simulator_.timeStepSize(); // Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES) - computeWellGroupThp(dt, deferred_logger); + bool well_group_thp_updated = computeWellGroupThp(dt, deferred_logger); constexpr int max_number_of_sub_iterations = 20; constexpr Scalar damping_factor = 0.1; + bool more_network_sub_update = false; for (int i = 0; i < max_number_of_sub_iterations; i++) { const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx, damping_factor); const Scalar network_imbalance = comm.max(local_network_imbalance); const auto& balance = this->schedule()[episodeIdx].network_balance(); constexpr Scalar relaxation_factor = 10.0; const Scalar tolerance = relax_network_tolerance ? relaxation_factor * balance.pressure_tolerance() : balance.pressure_tolerance(); - more_network_update = this->networkActive() && network_imbalance > tolerance; - if (!more_network_update) + more_network_sub_update = this->networkActive() && network_imbalance > tolerance; + if (!more_network_sub_update) break; } + more_network_update = more_network_sub_update || well_group_thp_updated; } bool changed_well_group = false; diff --git a/opm/simulators/wells/FractionCalculator.cpp b/opm/simulators/wells/FractionCalculator.cpp index 6af54ee58..2a7230c66 100644 --- a/opm/simulators/wells/FractionCalculator.cpp +++ b/opm/simulators/wells/FractionCalculator.cpp @@ -150,7 +150,7 @@ guideRate(const std::string& name, return WellGroupHelpers::getGuideRate(name, schedule_, well_state_, group_state_, report_step_, guide_rate_, target_, pu_); } else { - if ((groupControlledWells(name, always_included_child) > 0)) { + if (groupControlledWells(name, always_included_child) > 0) { if (is_producer_ && guide_rate_->has(name)) { return guide_rate_->get(name, target_, getGroupRateVector(name)); } else if (!is_producer_ && guide_rate_->has(name, injection_phase_)) { @@ -178,6 +178,7 @@ groupControlledWells(const std::string& group_name, well_state_, this->group_state_, this->summary_state_, + this->guide_rate_, report_step_, group_name, always_included_child, diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index b60ece73c..8407c0908 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -1058,7 +1058,6 @@ bruteForceBracketCommonTHP(const std::function& eq, for (int i = 0; i < sample_number + 1; ++i) { max_thp = min_thp + interval * i; eq_high = eq(max_thp); - // std::cout << "max_thp: " << max_thp/1E5 << " eq_high: " << eq_high << " eq_low: " << eq_low << std::endl; if (eq_high * eq_low <= 0.) { bracket_found = true; min_thp = max_thp - interval; diff --git a/opm/simulators/wells/WellGroupControls.cpp b/opm/simulators/wells/WellGroupControls.cpp index a91e482bc..5ed3f3d3d 100644 --- a/opm/simulators/wells/WellGroupControls.cpp +++ b/opm/simulators/wells/WellGroupControls.cpp @@ -529,17 +529,14 @@ getGroupProductionTargetRate(const Group& group, // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. const std::size_t num_ancestors = chain.size() - 1; Scalar target = orig_target; - std::cout << "target: " << target*86400 << " time: " << well_.currentStep() << " orig" << std::endl; for (std::size_t ii = 0; ii < num_ancestors; ++ii) { if ((ii == 0) || well_.guideRate()->has(chain[ii])) { // Apply local reductions only at the control level // (top) and for levels where we have a specified // group guide rate. target -= localReduction(chain[ii]); - std::cout << "ii: " << ii << " group: " << chain[ii] << " LocRed:" << localReduction(chain[ii])*86400 << " target: " << target*86400 << std::endl; } target *= localFraction(chain[ii+1]); - std::cout << "ii: " << ii << " group: " << chain[ii+1] << " LocFrac: " << localFraction(chain[ii+1]) << " target: " << target*86400 << std::endl; } // Avoid negative target rates coming from too large local reductions. const Scalar target_rate = std::max(Scalar(0.0), target / efficiencyFactor); @@ -636,29 +633,20 @@ getAutoChokeGroupProductionTargetRate(const std::string& name, if (!group.has_gpmaint_control(currentGroupControl)) ctrl = group.productionControls(summaryState); - // Scalar fr_true = fcalc.fraction("B1", "M5S", true); - // fr_true = fcalc.fraction("B1", "M5S", true); - // Scalar fr_false = fcalc.fraction("B1", "M5S", false); - // fr_false = fcalc.fraction("B1", "M5S", false); - // std::cout << "fr_true: " << fr_true << " fr_false: " << fr_false << std::endl; - const double orig_target = tcalc.groupTarget(ctrl, deferred_logger); const auto chain = WellGroupHelpers::groupChainTopBot(name, group.name(), schedule, reportStepIdx); // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. const std::size_t num_ancestors = chain.size() - 1; double target = orig_target; - std::cout << "target: " << target*86400 << " time: " << reportStepIdx << " modified" << std::endl; for (std::size_t ii = 0; ii < num_ancestors; ++ii) { if ((ii == 0) || guideRate->has(chain[ii])) { // Apply local reductions only at the control level // (top) and for levels where we have a specified // group guide rate. target -= localReduction(chain[ii]); - std::cout << "ii: " << ii << " group: " << chain[ii] << " LocRed:" << localReduction(chain[ii])*86400 << " target: " << target*86400 << std::endl; } target *= localFraction(chain[ii+1]); - std::cout << "ii: " << ii << " group: " << chain[ii+1] << " LocFrac: " << localFraction(chain[ii+1]) << " target: " << target*86400 << std::endl; } // Avoid negative target rates coming from too large local reductions. const double target_rate = std::max(0.0, target / efficiencyFactor); diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index bb6a133ff..755bd3057 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -404,7 +404,7 @@ updateGroupTargetReduction(const Group& group, const bool individual_control = (currentGroupControl != Group::InjectionCMode::FLD && currentGroupControl != Group::InjectionCMode::NONE); const int num_group_controlled_wells - = groupControlledWells(schedule, wellState, group_state, summaryState, reportStepIdx, subGroupName, "", !isInjector, phase); + = groupControlledWells(schedule, wellState, group_state, summaryState, &guide_rate, reportStepIdx, subGroupName, "", !isInjector, phase); if (individual_control || num_group_controlled_wells == 0) { groupTargetReduction[phase_pos] += subGroupEfficiency * sumWellSurfaceRates(subGroup, schedule, wellState, reportStepIdx, phase_pos, isInjector); @@ -420,8 +420,7 @@ updateGroupTargetReduction(const Group& group, const bool individual_control = (currentGroupControl != Group::ProductionCMode::FLD && currentGroupControl != Group::ProductionCMode::NONE); const int num_group_controlled_wells - = groupControlledWells(schedule, wellState, group_state, summaryState, reportStepIdx, subGroupName, "", !isInjector, /*injectionPhaseNotUsed*/Phase::OIL); - std::cout<< "Group: " << subGroupName << " NumGroupControlledWells: " << num_group_controlled_wells << std::endl; + = groupControlledWells(schedule, wellState, group_state, summaryState, &guide_rate, reportStepIdx, subGroupName, "", !isInjector, /*injectionPhaseNotUsed*/Phase::OIL); if (individual_control || num_group_controlled_wells == 0) { for (int phase = 0; phase < np; phase++) { groupTargetReduction[phase] @@ -436,7 +435,6 @@ updateGroupTargetReduction(const Group& group, } } } - std::cout << "A: Group: " << group.name() << ", currentGroupControl: " << Group::ProductionCMode2String(currentGroupControl) << ", GroupTargetReduction[1]: " << groupTargetReduction[1]*86400 << std::endl; } } @@ -476,9 +474,6 @@ updateGroupTargetReduction(const Group& group, groupTargetReduction[phase] -= ws.surface_rates[phase] * efficiency; } } - std::cout<< "B: " << " Group: " << group.name() << ", Well: " - << wellName << ", cmode: " << ws.production_cmode << ", GroupTargetReduction[1]: " - << groupTargetReduction[1]*86400 << std::endl; } } } @@ -1099,6 +1094,7 @@ groupControlledWells(const Schedule& schedule, const WellState& well_state, const GroupState& group_state, const SummaryState& summary_state, + const GuideRate* guideRate, const int report_step, const std::string& group_name, const std::string& always_included_child, @@ -1120,7 +1116,7 @@ groupControlledWells(const Schedule& schedule, if (included) { num_wells - += groupControlledWells(schedule, well_state, group_state, summary_state, report_step, child_group, always_included_child, is_production_group, injection_phase); + += groupControlledWells(schedule, well_state, group_state, summary_state, guideRate, report_step, child_group, always_included_child, is_production_group, injection_phase); } } for (const std::string& child_well : group.wells()) { @@ -1152,31 +1148,47 @@ groupControlledWells(const Schedule& schedule, const auto& control_group_name = control_group(group, group_state, report_step, schedule); const auto& control_group = schedule.getGroup(control_group_name, report_step); const auto& ctrl = control_group.productionControls(summary_state); - const auto& control_group_guide_rate = ctrl.guide_rate; const auto& control_group_cmode = ctrl.cmode; const auto& group_guide_rate = group.productionControls(summary_state).guide_rate; - Scalar gratTargetFromSales = 0.0; - if (group_state.has_grat_sales_target(control_group_name)) - gratTargetFromSales = group_state.grat_sales_target(control_group_name); + if (group_guide_rate > 0) { + // Guide rate is not default for the auto choke group + Scalar gratTargetFromSales = 0.0; + if (group_state.has_grat_sales_target(control_group_name)) + gratTargetFromSales = group_state.grat_sales_target(control_group_name); - std::vector resv_coeff(pu.num_phases, 1.0); - WGHelpers::TargetCalculator tcalc(control_group_cmode, - pu, - resv_coeff, - gratTargetFromSales, - group.name(), - group_state, - group.has_gpmaint_control(control_group_cmode)); - auto deferred_logger = Opm::DeferredLogger(); - const auto& control_group_target = tcalc.groupTarget(ctrl, deferred_logger); - // Target rate for the auto choke group - const Scalar target_rate = control_group_target * group_guide_rate / control_group_guide_rate; - const Scalar current_rate = tcalc.calcModeRateFromRates(rates); + std::vector resv_coeff(pu.num_phases, 1.0); + WGHelpers::TargetCalculator tcalc(control_group_cmode, + pu, + resv_coeff, + gratTargetFromSales, + group.name(), + group_state, + group.has_gpmaint_control(control_group_cmode)); + auto deferred_logger = Opm::DeferredLogger(); + const auto& control_group_target = tcalc.groupTarget(ctrl, deferred_logger); - if (current_rate < 0.99 * target_rate) - included = false; + // Calculates the guide rate of the parent group with control. + // It is allowed that the guide rate of this group is defaulted. The guide rate will be derived from the children groups + const auto& control_group_guide_rate = getGuideRate(control_group_name, + schedule, + well_state, + group_state, + report_step, + guideRate, + tcalc.guideTargetMode(), + pu); + + if (control_group_guide_rate > 0) { + // Target rate for the auto choke group + const Scalar target_rate = control_group_target * group_guide_rate / control_group_guide_rate; + const Scalar current_rate = tcalc.calcModeRateFromRates(rates); + + if (current_rate < target_rate) + included = false; + } + } } if (included) { @@ -1361,6 +1373,7 @@ checkGroupConstraintsProd(const std::string& name, wellState, group_state, summaryState, + guideRate, reportStepIdx, chain[ii], "", @@ -1381,7 +1394,6 @@ checkGroupConstraintsProd(const std::string& name, // we add a factor here to avoid switching due to numerical instability const Scalar factor = 1.01; if (currentRateFraction > (guiderateFraction * factor)) { - std::cout << guided_group << " localCurrentRate(guided_group): " << localCurrentRate(guided_group) << chain[ii-1] << " localCurrentRate(chain[ii-1]) " << localCurrentRate(chain[ii-1]) << std::endl; return std::make_pair(true, guiderateFraction / currentRateFraction); } } @@ -1540,6 +1552,7 @@ checkGroupConstraintsInj(const std::string& name, wellState, group_state, summaryState, + guideRate, reportStepIdx, chain[ii], "", diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 9687954b0..00db85713 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -249,6 +249,7 @@ public: const WellState& well_state, const GroupState& group_state, const SummaryState& summary_state, + const GuideRate* guideRate, const int report_step, const std::string& group_name, const std::string& always_included_child, diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index af95e9cd4..00983b534 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1482,7 +1482,6 @@ namespace Opm } else { ws.trivial_target = true; } - break; } case Well::ProducerCMode::CMODE_UNDEFINED: