From 30e1f5178c3e1aeb4dcaafbbd38bd624cb91aa60 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 3 Oct 2023 13:28:23 +0200 Subject: [PATCH 1/8] Automatic choke --- opm/simulators/wells/BlackoilWellModel.hpp | 2 + .../wells/BlackoilWellModelGeneric.cpp | 4 +- .../wells/BlackoilWellModel_impl.hpp | 85 ++++++++++++++++++- opm/simulators/wells/GroupState.cpp | 17 ++++ opm/simulators/wells/GroupState.hpp | 5 ++ opm/simulators/wells/WellAssemble.cpp | 36 +++++--- opm/simulators/wells/WellBhpThpCalculator.hpp | 5 ++ opm/simulators/wells/WellGroupHelpers.cpp | 11 ++- opm/simulators/wells/WellInterface_impl.hpp | 41 +++++---- 9 files changed, 175 insertions(+), 31 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 0fe1ef121..21496b72e 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -465,6 +465,8 @@ template class WellContributions; const double dt, DeferredLogger& local_deferredLogger); + void computeWellGroupThp(DeferredLogger& local_deferredLogger); + /// Update rank's notion of intersecting wells and their /// associate solution variables. /// diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 32b9ea5ef..0d60fd97b 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1363,7 +1363,6 @@ updateNetworkPressures(const int reportStepIdx) } 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. @@ -1380,6 +1379,9 @@ updateNetworkPressures(const int reportStepIdx) 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 6354bfcf4..2af37e1b6 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -28,6 +28,7 @@ #endif #include +#include #include #include @@ -40,6 +41,9 @@ #include #include #include +#include +#include +#include #include #include @@ -1162,6 +1166,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; @@ -1243,7 +1250,82 @@ namespace Opm { return {more_network_update, well_group_control_changed}; } + //PJPE: This function is to be used for well groups in an extended network that act as a subsea manifold + // 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 + // subsea manifolds + template + void + BlackoilWellModel:: + computeWellGroupThp(DeferredLogger& local_deferredLogger) + { + const int reportStepIdx = this->ebosSimulator_.episodeIndex(); + const auto& network = schedule()[reportStepIdx].network(); + if (!network.active()) { + return; + } + + auto& well_state = this->wellState(); + auto& group_state = this->groupState(); + + for (const std::string& nodeName : network.node_names()) { + const bool has_choke = network.node(nodeName).as_choke(); + if (has_choke) { + const auto& summary_state = this->ebosSimulator_.vanguard().summaryState(); + const Group& group = schedule().getGroup(nodeName, reportStepIdx); + std::optional ctrl = group.productionControls(summary_state); + const auto cmode = ctrl->cmode; + const auto pu = this->phase_usage_; + //PJPE: conversion factor res<-> surface rates TODO: to be calculated + std::vector resv_coeff(pu.num_phases, 1.0); + double gratTargetFromSales = 0.0; //PJPE: check this + WellGroupHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff, + gratTargetFromSales, nodeName, group_state, + group.has_gpmaint_control(cmode)); + double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); + //PJPE: TODO: include efficiency factor + auto mismatch = [&] (auto well_group_thp) { + double well_group_rate(0.0); + for (auto& well : this->well_container_) { + std::string well_name = well->name(); + if (group.hasWell(well_name)) { + well->setDynamicThpLimit(well_group_thp); + well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); + auto& ws = well_state.well(well_name); + well_group_rate += -tcalc.calcModeRateFromRates(ws.surface_rates); + } + } + return well_group_rate - orig_target; + }; + + // PJPE: TODO what settings to take here? + const double tolerance = 1E-4; + const int max_iteration_solve = 100; + int iteration = 0; + const std::array range {20E5, 40E5}; + double low(range[0]); + double high(range[1]); + + // PJPE: TODO use bisectBracket first and use bruteForceBracket as fallback option. + // see computeBhpAtThpLimit() + WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger); + + // double aa = mismatch(low); + // double bb = mismatch(high); + // std::cout << " low = " << aa << ", high = " << bb << std::endl; + + double well_group_thp = RegulaFalsiBisection:: + solve(mismatch, low, high, max_iteration_solve, tolerance, iteration); + + double check = mismatch(well_group_thp); + std::cout << " mismatch = " << check << ", well_group_thp = " << well_group_thp << std::endl; + + //PJPE: use the common THP in computeNetworkPressures + group_state.update_well_group_thp(nodeName, well_group_thp); + } + } + } template void @@ -2374,9 +2456,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 85167096a..a46ed4ef1 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -101,6 +101,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 bool GroupState:: diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index 15f90dab6..97ab1d0f0 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -49,6 +49,9 @@ public: const std::vector& rates); const std::vector& production_rates(const std::string& gname) const; + void update_well_group_thp(const std::string& gname, const double& thp); + double well_group_thp(const std::string& gname) const; + bool has_production_reduction_rates(const std::string& gname) const; void update_production_reduction_rates(const std::string& gname, const std::vector& rates); @@ -175,6 +178,7 @@ public: serializer(num_phases); serializer(m_production_rates); serializer(production_controls); + serializer(group_thp); serializer(prod_red_rates); serializer(inj_red_rates); serializer(inj_surface_rates); @@ -199,6 +203,7 @@ private: std::map inj_vrep_rate; std::map m_grat_sales_target; std::map m_gpmaint_target; + std::map group_thp; std::map, Group::InjectionCMode> injection_controls; WellContainer gpmaint_state; diff --git a/opm/simulators/wells/WellAssemble.cpp b/opm/simulators/wells/WellAssemble.cpp index e2994a0e6..83e77c3cc 100644 --- a/opm/simulators/wells/WellAssemble.cpp +++ b/opm/simulators/wells/WellAssemble.cpp @@ -35,6 +35,7 @@ #include #include #include +#include #include #include @@ -162,17 +163,30 @@ assembleControlEqProd(const WellState& well_state, well_.rateConverter().calcCoeff(id, region, coeff); }; - WellGroupControls(well_).getGroupProductionControl(group, - well_state, - group_state, - schedule, - summaryState, - bhp, - active_rates, - rCoeff, - efficiencyFactor, - control_eq, - deferred_logger); + + bool has_choke(false); + const std::size_t report_step = well_.currentStep(); + auto& network = schedule[report_step].network(); + if (network.active()) + has_choke = network.node(group.name()).as_choke(); + + if(!has_choke) { + WellGroupControls(well_).getGroupProductionControl(group, + well_state, + group_state, + schedule, + summaryState, + bhp, + active_rates, + rCoeff, + efficiencyFactor, + control_eq, + deferred_logger); + + } else { + //PJPE: the group is a subsea manifold: wells are operated on a common THP + control_eq = bhp - bhp_from_thp(); + } break; } case Well::ProducerCMode::CMODE_UNDEFINED: { diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index 10c059e01..f0c80e864 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -118,6 +118,11 @@ public: getFloIPR(const WellState& well_state, const Well& well, const SummaryState& summary_state) const; + //! \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); private: //! \brief Compute BHP from THP limit for an injector - implementation. diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 472a80d63..5bf14af7b 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -41,6 +41,7 @@ #include #include #include +#include #include #include @@ -935,7 +936,15 @@ computeNetworkPressures(const Network::ExtNetwork& network, #endif } else { // Table number specified as 9999 in the deck, no pressure loss. - node_pressures[node] = up_press; + if (network.node(node).as_choke()){ + // Node pressure is set to the common THP of the wells. + // The choke pressure must be non-negative therefore the node pressure of + // the auto-choke node must be larger or equal to the pressure of the uptree node of its branch + const auto group_thp = group_state.well_group_thp(node); + node_pressures[node] = group_thp >= up_press ? group_thp : up_press; + } else { + node_pressures[node] = up_press; + } } } } diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index ec451bc56..9818bf162 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1401,25 +1401,36 @@ namespace Opm assert(well.isAvailableForGroupControl()); const auto& group = schedule.getGroup(well.groupName(), this->currentStep()); const Scalar efficiencyFactor = well.getEfficiencyFactor(); - Scalar scale = this->getGroupProductionTargetRate(group, - well_state, - group_state, - schedule, - summaryState, - efficiencyFactor, - deferred_logger); - // we don't want to scale with zero and get zero rates. - if (scale > 0) { - for (int p = 0; pcurrentStep(); + auto& network = schedule[report_step].network(); + if (network.active()) + has_choke = network.node(group.name()).as_choke(); + if (!has_choke) { + Scalar scale = this->getGroupProductionTargetRate(group, + well_state, + group_state, + schedule, + summaryState, + efficiencyFactor, + deferred_logger); + + // we don't want to scale with zero and get zero rates. + if (scale > 0) { + for (int p = 0; p Date: Thu, 2 Nov 2023 10:59:16 +0100 Subject: [PATCH 2/8] autochoke producers control mode set to THP --- .../wells/BlackoilWellModel_impl.hpp | 55 ++++++++++++++----- opm/simulators/wells/WellAssemble.cpp | 34 ++++-------- opm/simulators/wells/WellInterface_impl.hpp | 1 + 3 files changed, 53 insertions(+), 37 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 2af37e1b6..dabb265ae 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1253,7 +1253,6 @@ namespace Opm { //PJPE: This function is to be used for well groups in an extended network that act as a subsea manifold // 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 - // subsea manifolds template void BlackoilWellModel:: @@ -1299,27 +1298,55 @@ namespace Opm { return well_group_rate - orig_target; }; + const std::array range {20E5, 40E5}; + double low, high; + + // // // trying to use bisect way to locate a bracket + // std::optional approximate_solution; + // bool finding_bracket = WellBhpThpCalculator::bisectBracketModified(mismatch, range, low, high, approximate_solution, local_deferredLogger); + // if (approximate_solution.has_value()) + // std::cout << "Approximate solution: " << *approximate_solution << std::endl; + + bool finding_bracket = false; + if (!finding_bracket) { + local_deferredLogger.debug(" Trying the brute force search to bracket the bhp for last attempt "); + finding_bracket = WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger); + } + + if (!finding_bracket) { + local_deferredLogger.warning("Bracketing THP calculation failed"); + // return std::nullopt; + } + // PJPE: TODO what settings to take here? const double tolerance = 1E-4; const int max_iteration_solve = 100; int iteration = 0; - const std::array range {20E5, 40E5}; - double low(range[0]); - double high(range[1]); - - // PJPE: TODO use bisectBracket first and use bruteForceBracket as fallback option. - // see computeBhpAtThpLimit() - WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger); - - // double aa = mismatch(low); - // double bb = mismatch(high); - // std::cout << " low = " << aa << ", high = " << bb << std::endl; - double well_group_thp = RegulaFalsiBisection:: solve(mismatch, low, high, max_iteration_solve, tolerance, iteration); + // std::cout << "Bracket = [" << low << ", " << high <<"]" << std::endl; + // if (low < 2.5E6){ + // double pr(0); + // for (int i = 0; i < 1000; ++i){ + // pr = (i/10.0 + 1) * 1.0E5 ; + // // std::cout << " thp = " << pr[i] << ", mismatch = " << mismatch(pr[i]) << std::endl; + // std::cout << pr << ", " << mismatch(pr) << std::endl; + // } + // return; + // } + + for (auto& well : this->well_container_) { + std::string well_name = well->name(); + if (group.hasWell(well_name)) { + well->setDynamicThpLimit(well_group_thp); + auto& ws = well_state.well(well_name); + ws.thp = well_group_thp; + } + } + double check = mismatch(well_group_thp); - std::cout << " mismatch = " << check << ", well_group_thp = " << well_group_thp << std::endl; + std::cout << "mismatch = " << check << ", well_group_thp = " << well_group_thp << ", iteration = " << iteration << std::endl; //PJPE: use the common THP in computeNetworkPressures group_state.update_well_group_thp(nodeName, well_group_thp); diff --git a/opm/simulators/wells/WellAssemble.cpp b/opm/simulators/wells/WellAssemble.cpp index 83e77c3cc..46314a8fd 100644 --- a/opm/simulators/wells/WellAssemble.cpp +++ b/opm/simulators/wells/WellAssemble.cpp @@ -163,30 +163,18 @@ assembleControlEqProd(const WellState& well_state, well_.rateConverter().calcCoeff(id, region, coeff); }; - - bool has_choke(false); - const std::size_t report_step = well_.currentStep(); - auto& network = schedule[report_step].network(); - if (network.active()) - has_choke = network.node(group.name()).as_choke(); - if(!has_choke) { - WellGroupControls(well_).getGroupProductionControl(group, - well_state, - group_state, - schedule, - summaryState, - bhp, - active_rates, - rCoeff, - efficiencyFactor, - control_eq, - deferred_logger); - - } else { - //PJPE: the group is a subsea manifold: wells are operated on a common THP - control_eq = bhp - bhp_from_thp(); - } + WellGroupControls(well_).getGroupProductionControl(group, + well_state, + group_state, + schedule, + summaryState, + bhp, + active_rates, + rCoeff, + efficiencyFactor, + control_eq, + deferred_logger); break; } case Well::ProducerCMode::CMODE_UNDEFINED: { diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 9818bf162..0d447b6e6 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1428,6 +1428,7 @@ namespace Opm } else { // PJPE: the group is a subsea manifold.Guide rates to be ignored. // The wells of the group are to be operated on a common THP (= manifold node pressure) + ws.production_cmode = Well::ProducerCMode::THP; } break; } From 36dcf233c533be02991bfa60223990ede4a71f56 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 7 Nov 2023 22:24:24 +0100 Subject: [PATCH 3/8] simplyfying code as some logic is moved to the parser --- .../wells/BlackoilWellModelGeneric.cpp | 3 -- .../wells/BlackoilWellModel_impl.hpp | 24 +++++------ opm/simulators/wells/WellInterface_impl.hpp | 40 +++++++------------ 3 files changed, 23 insertions(+), 44 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 0d60fd97b..6f6eb88d3 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1379,9 +1379,6 @@ updateNetworkPressures(const int reportStepIdx) 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 dabb265ae..3665ed638 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1298,24 +1298,15 @@ namespace Opm { return well_group_rate - orig_target; }; - const std::array range {20E5, 40E5}; + const std::array range {1E5, 150E5}; //PJPE what lower/upper bound to be taken here? double low, high; - // // // trying to use bisect way to locate a bracket - // std::optional approximate_solution; - // bool finding_bracket = WellBhpThpCalculator::bisectBracketModified(mismatch, range, low, high, approximate_solution, local_deferredLogger); - // if (approximate_solution.has_value()) - // std::cout << "Approximate solution: " << *approximate_solution << std::endl; - - bool finding_bracket = false; - if (!finding_bracket) { - local_deferredLogger.debug(" Trying the brute force search to bracket the bhp for last attempt "); - finding_bracket = WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger); - } - - if (!finding_bracket) { + local_deferredLogger.debug(" Trying the brute force search to bracket the bhp for last attempt "); + const bool finding_bracket = WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger); + if (finding_bracket) { + std::cout << "low: " << low << " high: " << high << std::endl; + } else { local_deferredLogger.warning("Bracketing THP calculation failed"); - // return std::nullopt; } // PJPE: TODO what settings to take here? @@ -2482,6 +2473,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 0d447b6e6..098c6ff1c 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1401,35 +1401,23 @@ namespace Opm assert(well.isAvailableForGroupControl()); const auto& group = schedule.getGroup(well.groupName(), this->currentStep()); const Scalar efficiencyFactor = well.getEfficiencyFactor(); + Scalar scale = this->getGroupProductionTargetRate(group, + well_state, + group_state, + schedule, + summaryState, + efficiencyFactor, + deferred_logger); - bool has_choke(false); - const std::size_t report_step = this->currentStep(); - auto& network = schedule[report_step].network(); - if (network.active()) - has_choke = network.node(group.name()).as_choke(); - if (!has_choke) { - Scalar scale = this->getGroupProductionTargetRate(group, - well_state, - group_state, - schedule, - summaryState, - 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 4/8] allow individual well constraints --- opm/simulators/wells/BlackoilWellModel.hpp | 2 +- .../wells/BlackoilWellModel_impl.hpp | 104 +++++++++++++----- opm/simulators/wells/WellBhpThpCalculator.cpp | 5 + opm/simulators/wells/WellBhpThpCalculator.hpp | 6 + opm/simulators/wells/WellGroupHelpers.cpp | 2 +- opm/simulators/wells/WellInterface.hpp | 9 ++ 6 files changed, 99 insertions(+), 29 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 21496b72e..bbc3092a3 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -465,7 +465,7 @@ template class WellContributions; const double dt, DeferredLogger& local_deferredLogger); - void computeWellGroupThp(DeferredLogger& local_deferredLogger); + void 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 3665ed638..ad869fd9c 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1167,8 +1167,7 @@ namespace Opm { 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); - + computeWellGroupThp(dt,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; @@ -1181,6 +1180,7 @@ namespace Opm { if (this->terminal_output_ && (network_update_iteration == iteration_to_relax) ) { local_deferredLogger.info(" we begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) + " iterations "); } + const bool relax_network_balance = network_update_iteration >= iteration_to_relax; std::tie(do_network_update, well_group_control_changed) = updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, dt,local_deferredLogger); @@ -1256,7 +1256,7 @@ namespace Opm { template void BlackoilWellModel:: - computeWellGroupThp(DeferredLogger& local_deferredLogger) + computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger) { const int reportStepIdx = this->ebosSimulator_.episodeIndex(); const auto& network = schedule()[reportStepIdx].network(); @@ -1284,61 +1284,111 @@ namespace Opm { group.has_gpmaint_control(cmode)); double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); //PJPE: TODO: include efficiency factor + auto mismatch = [&] (auto well_group_thp) { + // auto well_state_copy = this->wellState(); double well_group_rate(0.0); for (auto& well : this->well_container_) { std::string well_name = well->name(); if (group.hasWell(well_name)) { well->setDynamicThpLimit(well_group_thp); - well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); + const Well& well_ecl = wells_ecl_[well->indexOfWell()]; + const auto inj_controls = Well::InjectionControls(0); + const auto prod_controls = well_ecl.productionControls(summary_state); + well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger); + // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); auto& ws = well_state.well(well_name); - well_group_rate += -tcalc.calcModeRateFromRates(ws.surface_rates); + double rate = -tcalc.calcModeRateFromRates(ws.surface_rates); + std::cout << "well: " << well_name << "rate: " << rate << std::endl; + well_group_rate += rate; } } + // this->wellState() = well_state; return well_group_rate - orig_target; }; - const std::array range {1E5, 150E5}; //PJPE what lower/upper bound to be taken here? + const auto upbranch = network.uptree_branch(nodeName); + const auto it = node_pressures_.find((*upbranch).uptree_node()); + const double nodal_pressure = it->second; + const std::array range {1.0E5, 150.0E5}; //PJPE what lower/upper bound to be taken here? double low, high; + low = 3527000.0; + high = 3576666.7; + std::cout << "low: " << low << " high: " << high << std::endl; + double low_v = mismatch(low); + std::cout << "1: low_value: " << low_v << std::endl; + low_v = mismatch(low); + high_v = mismatch(high); + std::cout << "2: low_value: " << low_v << " high_value: " << high_v << std::endl; - local_deferredLogger.debug(" Trying the brute force search to bracket the bhp for last attempt "); + + local_deferredLogger.debug(" Trying the brute force search to bracket the common THP "); const bool finding_bracket = WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger); - if (finding_bracket) { - std::cout << "low: " << low << " high: " << high << std::endl; + std::cout << "low: " << low << " high: " << high << std::endl; + double low_value = mismatch(low); + double high_value = mismatch(high); + std::cout << "low_value: " << mismatch(low) << " high_value: " << mismatch(high) << std::endl; + + double well_group_thp = 0.0; + if ((low_value * high_value >= 0.0) && (low >= nodal_pressure) ) { + // In case no bracket is found + const double limit = 0.0003; // PJPE what limit to be used? + double abs_low = std::fabs(low_value); + double abs_high = std::fabs(high_value); + if (std::min(abs_low, abs_high) < limit) { + well_group_thp = abs_low < abs_high ? low : high; + local_deferredLogger.warning("Common THP not solved precisely"); + } else { + local_deferredLogger.warning("Common THP solve failed due to bracketing failure"); + well_group_thp = nodal_pressure; + } + } else if (high <= nodal_pressure) { + local_deferredLogger.warning("Common THP set to nodal pressure"); + well_group_thp = nodal_pressure; } else { - local_deferredLogger.warning("Bracketing THP calculation failed"); + // PJPE: TODO what settings to take here? + const double tolerance = 1E-4; + const int max_iteration_solve = 20; + int iteration = 0; + well_group_thp = RegulaFalsiBisection:: + solve(mismatch, low, high, max_iteration_solve, tolerance, iteration); + std::cout << "Bracket = [" << low << ", " << high <<"] " << ", iteration = " << iteration << std::endl; + double check = mismatch(well_group_thp); + std::cout << "mismatch = " << check << ", well_group_thp = " << well_group_thp << std::endl; } - // PJPE: TODO what settings to take here? - const double tolerance = 1E-4; - const int max_iteration_solve = 100; - int iteration = 0; - double well_group_thp = RegulaFalsiBisection:: - solve(mismatch, low, high, max_iteration_solve, tolerance, iteration); + // well_group_thp = std::max(well_group_thp, nodal_pressure); + std::cout << "well_group_thp = " << well_group_thp << ", nodal_pressure = " << nodal_pressure << std::endl; - // std::cout << "Bracket = [" << low << ", " << high <<"]" << std::endl; - // if (low < 2.5E6){ + // if (well_group_thp < 2.5E6){ + // double min = 24.0E5; + // double max = 34.0E5; // double pr(0); // for (int i = 0; i < 1000; ++i){ - // pr = (i/10.0 + 1) * 1.0E5 ; + // // pr = (i/10.0 + 1) * 1.0E5 ; + // pr = min + i * (max-min)/1000.0; // // std::cout << " thp = " << pr[i] << ", mismatch = " << mismatch(pr[i]) << std::endl; // std::cout << pr << ", " << mismatch(pr) << std::endl; // } // return; // } - for (auto& well : this->well_container_) { - std::string well_name = well->name(); - if (group.hasWell(well_name)) { + std::string well_name = well->name(); + if (group.hasWell(well_name)) { + auto& ws = well_state.well(well_name); + std::cout << "Cmode of Well: " << well_name << " is " << ws.production_cmode <setDynamicThpLimit(well_group_thp); - auto& ws = well_state.well(well_name); ws.thp = well_group_thp; + const Well& well_ecl = wells_ecl_[well->indexOfWell()]; + const auto inj_controls = Well::InjectionControls(0); + const auto prod_controls = well_ecl.productionControls(summary_state); + // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); + well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger); } } - - double check = mismatch(well_group_thp); - std::cout << "mismatch = " << check << ", well_group_thp = " << well_group_thp << ", iteration = " << iteration << std::endl; - + } + //PJPE: use the common THP in computeNetworkPressures group_state.update_well_group_thp(nodeName, well_group_thp); } diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 4b3ccadc2..bb6fcd9ba 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -889,6 +889,11 @@ bruteForceBracket(const std::function& eq, eq_low = eq_high; } if (bracket_found) { + // std::cout << "BFB: low: " << low << " high: " << high << std::endl; + // std::cout << "BFB: low_value: " << eq_low << " high_value: " << eq_high << std::endl; + // double low_value = eq(low); + // double high_value = eq(high); + // std::cout << "CHK: low_value: " << low_value<< " high_value: " << high_value << std::endl; deferred_logger.debug( " brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) + " high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high)); diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index f0c80e864..24806d31e 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -124,6 +124,12 @@ public: double& low, double& high, DeferredLogger& deferred_logger); + //! \brief Find limits using brute-force solver. + static bool bruteForceBracketHigh(const std::function& eq, + const std::array& range, + double& low, double& high, + DeferredLogger& deferred_logger); + private: //! \brief Compute BHP from THP limit for an injector - implementation. template diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 5bf14af7b..6ca7bd441 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -943,7 +943,7 @@ computeNetworkPressures(const Network::ExtNetwork& network, const auto group_thp = group_state.well_group_thp(node); node_pressures[node] = group_thp >= up_press ? group_thp : up_press; } else { - node_pressures[node] = up_press; + node_pressures[node] = up_press; } } } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 957bbfe6f..bd78bdcfd 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -374,6 +374,15 @@ public: void updateConnectionTransmissibilityFactor(const Simulator& simulator, SingleWellState& ws) const; + SingleWellState& ws) const; + + virtual bool iterateWellEqWithSwitching(const Simulator& simulator, + const double dt, + const WellInjectionControls& inj_controls, + const WellProductionControls& prod_controls, + WellState& well_state, + const GroupState& group_state, + DeferredLogger& deferred_logger) = 0; protected: // simulation parameters const ModelParameters& param_; From 1ddf675cfd8ac0aec29e74125cda847f45bd7eec Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 22 Jan 2024 13:47:06 +0100 Subject: [PATCH 5/8] moved common thp calculation to updateWellControls --- opm/simulators/wells/BlackoilWellModel.hpp | 1 + .../wells/BlackoilWellModel_impl.hpp | 172 +++++++++++------- opm/simulators/wells/WellBhpThpCalculator.cpp | 52 +++++- opm/simulators/wells/WellBhpThpCalculator.hpp | 25 +-- opm/simulators/wells/WellGroupHelpers.cpp | 2 +- opm/simulators/wells/WellInterface.hpp | 13 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 7 + opm/simulators/wells/WellInterfaceGeneric.hpp | 1 + 8 files changed, 185 insertions(+), 88 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index bbc3092a3..b65004ea8 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -413,6 +413,7 @@ template class WellContributions; Scalar gravity_{}; std::vector depth_{}; bool alternative_well_rate_init_{}; + double well_group_thp_prev_{0.0}; std::unique_ptr rateConverter_{}; std::map> regionalAveragePressureCalculator_{}; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index ad869fd9c..2468dc96a 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -49,9 +49,12 @@ #include #include +<<<<<<< HEAD #if COMPILE_BDA_BRIDGE #include #endif +======= +>>>>>>> f92db77d3 (moved common thp calculation to updateWellControls) #if HAVE_MPI #include @@ -1166,8 +1169,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(dt,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; @@ -1250,7 +1251,7 @@ namespace Opm { return {more_network_update, well_group_control_changed}; } - //PJPE: This function is to be used for well groups in an extended network that act as a subsea manifold + // This function is to be used for well groups in an extended network that act as a subsea manifold // 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 @@ -1276,8 +1277,9 @@ namespace Opm { std::optional ctrl = group.productionControls(summary_state); const auto cmode = ctrl->cmode; const auto pu = this->phase_usage_; - //PJPE: conversion factor res<-> surface rates TODO: to be calculated + //PJPE: conversion factor res<-> surface rates TODO: to be calculated std::vector resv_coeff(pu.num_phases, 1.0); + // rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS. double gratTargetFromSales = 0.0; //PJPE: check this WellGroupHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff, gratTargetFromSales, nodeName, group_state, @@ -1286,92 +1288,116 @@ namespace Opm { //PJPE: TODO: include efficiency factor auto mismatch = [&] (auto well_group_thp) { - // auto well_state_copy = this->wellState(); double well_group_rate(0.0); + double rate(0.0); for (auto& well : this->well_container_) { std::string well_name = well->name(); + auto& ws = well_state.well(well_name); if (group.hasWell(well_name)) { + rate = -tcalc.calcModeRateFromRates(ws.surface_rates); + std::cout << "(1) well: " << well_name << ", rate: " << rate << ", thp: " << ws.thp << ", Cmode: " << ws.production_cmode << std::endl; well->setDynamicThpLimit(well_group_thp); + const Well& well_ecl = wells_ecl_[well->indexOfWell()]; const auto inj_controls = Well::InjectionControls(0); const auto prod_controls = well_ecl.productionControls(summary_state); - well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger); - // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); - auto& ws = well_state.well(well_name); - double rate = -tcalc.calcModeRateFromRates(ws.surface_rates); - std::cout << "well: " << well_name << "rate: " << rate << std::endl; + well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false); + // May be used instead if the wells of the sub-sea manifold has only THP constraints + // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); + + rate = -tcalc.calcModeRateFromRates(ws.surface_rates); + if (ws.production_cmode == Opm::WellProducerCMode::THP) + ws.thp = well_group_thp; + + std::cout << "(2) well: " << well_name << ", rate: " << rate << ", thp: " << ws.thp << ", Cmode: " << ws.production_cmode << std::endl; well_group_rate += rate; } } - // this->wellState() = well_state; - return well_group_rate - orig_target; + double diff = well_group_rate - orig_target; // (well_group_rate - orig_target)/orig_target; + std::cout << "well_group_thp: " << well_group_thp << " mismatch: " << diff << "\n" << std::endl; + return diff; }; + double well_group_rate(0.0); + double rate(0.0); + double thp(0.0); + double min_thp(1.0E8); + double max_thp(0.0); + for (auto& well : this->well_container_) { + std::string well_name = well->name(); + if (group.hasWell(well_name)) { + auto& ws = well_state.well(well_name); + rate = -tcalc.calcModeRateFromRates(ws.surface_rates); + well_group_rate += rate; + thp = ws.thp; + min_thp = std::min(min_thp, thp); + max_thp = std::max(max_thp, thp); + } + } + const auto upbranch = network.uptree_branch(nodeName); const auto it = node_pressures_.find((*upbranch).uptree_node()); const double nodal_pressure = it->second; - const std::array range {1.0E5, 150.0E5}; //PJPE what lower/upper bound to be taken here? - double low, high; - low = 3527000.0; - high = 3576666.7; - std::cout << "low: " << low << " high: " << high << std::endl; - double low_v = mismatch(low); - std::cout << "1: low_value: " << low_v << std::endl; - low_v = mismatch(low); - high_v = mismatch(high); - std::cout << "2: low_value: " << low_v << " high_value: " << high_v << std::endl; + double well_group_thp = nodal_pressure; + if ((reportStepIdx <= 1) || (well_group_rate > 0.9*orig_target)) { + std::array range; // = {1.0E5, 150.0E5}; //PJPE what lower/upper bound to be taken here? + if (well_group_thp_prev_ == 0) { + range = {0.9*min_thp, 1.1*max_thp}; + } + else { + range = {0.9* this->well_group_thp_prev_, 1.1*this->well_group_thp_prev_}; + } - local_deferredLogger.debug(" Trying the brute force search to bracket the common THP "); - const bool finding_bracket = WellBhpThpCalculator::bruteForceBracket(mismatch, range, low, high, local_deferredLogger); - std::cout << "low: " << low << " high: " << high << std::endl; - double low_value = mismatch(low); - double high_value = mismatch(high); - std::cout << "low_value: " << mismatch(low) << " high_value: " << mismatch(high) << std::endl; + double low, high; + // For testing + // low = 3.42767e+06;// 3.32833e+06; + // // high = 3.42767e+06;// 3.32833e+06; + // // std::cout << "low: " << low << " high: " << high << std::endl; + // double low_value_test = mismatch(low); + // low_value_test = mismatch(low); + // low_value_test = mismatch(low); + // std::cout << "low_value: " << low_value_test << " high_value: " << high_value_test << std::endl; - double well_group_thp = 0.0; - if ((low_value * high_value >= 0.0) && (low >= nodal_pressure) ) { - // In case no bracket is found - const double limit = 0.0003; // PJPE what limit to be used? - double abs_low = std::fabs(low_value); - double abs_high = std::fabs(high_value); - if (std::min(abs_low, abs_high) < limit) { - well_group_thp = abs_low < abs_high ? low : high; - local_deferredLogger.warning("Common THP not solved precisely"); + std::optional approximate_solution; + const double limit = 1.0E-3;// 1.0E-4; + local_deferredLogger.debug(" Using brute force search to bracket the common THP "); + std::cout << " Using brute force search to bracket the common THP " << std::endl; + const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, limit, local_deferredLogger); + + double low_value = mismatch(low); + double high_value = mismatch(high); + + if (approximate_solution.has_value()) { + well_group_thp = *approximate_solution; + double check = mismatch(well_group_thp); + local_deferredLogger.debug("Approximate value found: mismatch = " + std::to_string(check) + ", well_group_thp = " + std::to_string(well_group_thp)); + } else if (finding_bracket) { + // PJPE: TODO what settings to take here? + std::cout << "Found bracket: [" << low << ", " << high << "]; mismatch: [" << low_value << ", "<< high_value << "]" << std::endl; + const double tolerance = 1.0E-3; + const int max_iteration_solve = 100; + int iteration = 0; + well_group_thp = RegulaFalsiBisection:: + solve(mismatch, low, high, max_iteration_solve, tolerance, iteration); + local_deferredLogger.debug( " bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " + + "mismatch = [" + std::to_string(low_value) + ", " + std::to_string(high_value) + "], " + + "iteration = " + std::to_string(iteration)); + double check = mismatch(well_group_thp); + local_deferredLogger.debug(" mismatch = " + std::to_string(check) + ", well_group_thp = " + std::to_string(well_group_thp)); } else { local_deferredLogger.warning("Common THP solve failed due to bracketing failure"); well_group_thp = nodal_pressure; } - } else if (high <= nodal_pressure) { - local_deferredLogger.warning("Common THP set to nodal pressure"); - well_group_thp = nodal_pressure; + + std::cout << "well_group_thp = " << well_group_thp << ", nodal_pressure = " << nodal_pressure << std::endl; + } else { - // PJPE: TODO what settings to take here? - const double tolerance = 1E-4; - const int max_iteration_solve = 20; - int iteration = 0; - well_group_thp = RegulaFalsiBisection:: - solve(mismatch, low, high, max_iteration_solve, tolerance, iteration); - std::cout << "Bracket = [" << low << ", " << high <<"] " << ", iteration = " << iteration << std::endl; - double check = mismatch(well_group_thp); - std::cout << "mismatch = " << check << ", well_group_thp = " << well_group_thp << std::endl; + well_group_thp = nodal_pressure; } - // well_group_thp = std::max(well_group_thp, nodal_pressure); - std::cout << "well_group_thp = " << well_group_thp << ", nodal_pressure = " << nodal_pressure << std::endl; + well_group_thp = std::max(well_group_thp, nodal_pressure); - // if (well_group_thp < 2.5E6){ - // double min = 24.0E5; - // double max = 34.0E5; - // double pr(0); - // for (int i = 0; i < 1000; ++i){ - // // pr = (i/10.0 + 1) * 1.0E5 ; - // pr = min + i * (max-min)/1000.0; - // // std::cout << " thp = " << pr[i] << ", mismatch = " << mismatch(pr[i]) << std::endl; - // std::cout << pr << ", " << mismatch(pr) << std::endl; - // } - // return; - // } for (auto& well : this->well_container_) { std::string well_name = well->name(); if (group.hasWell(well_name)) { @@ -1380,16 +1406,19 @@ namespace Opm { if (ws.production_cmode == Opm::WellProducerCMode::THP){ well->setDynamicThpLimit(well_group_thp); ws.thp = well_group_thp; - const Well& well_ecl = wells_ecl_[well->indexOfWell()]; - const auto inj_controls = Well::InjectionControls(0); - const auto prod_controls = well_ecl.productionControls(summary_state); - // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); - well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger); } + // const Well& well_ecl = wells_ecl_[well->indexOfWell()]; + // const auto inj_controls = Well::InjectionControls(0); + // const auto prod_controls = well_ecl.productionControls(summary_state); + // // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); + // well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state_copy, group_state, local_deferredLogger); + if (ws.production_cmode == Opm::WellProducerCMode::THP) + ws.thp = well_group_thp; } } - - //PJPE: use the common THP in computeNetworkPressures + well_group_thp_prev_ = well_group_thp; + + // Use the common THP in computeNetworkPressures group_state.update_well_group_thp(nodeName, well_group_thp); } } @@ -2055,6 +2084,9 @@ namespace Opm { // network related bool more_network_update = false; if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) { + const double dt = this->ebosSimulator_.timeStepSize(); + // Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES) + computeWellGroupThp(dt, deferred_logger); const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx); const Scalar network_imbalance = comm.max(local_network_imbalance); const auto& balance = this->schedule()[episodeIdx].network_balance(); diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index bb6fcd9ba..9017bedb2 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -875,9 +875,15 @@ bruteForceBracket(const std::function& eq, low = range[0]; high = range[1]; const int sample_number = 200; +<<<<<<< HEAD const Scalar interval = (high - low) / sample_number; Scalar eq_low = eq(low); Scalar eq_high = 0.0; +======= + const double interval = (high - low) / sample_number; + double eq_low = eq(low); + double eq_high = 0.0; +>>>>>>> 3ce8d7eec (moved common thp calculation to updateWellControls) for (int i = 0; i < sample_number + 1; ++i) { high = range[0] + interval * i; eq_high = eq(high); @@ -889,11 +895,6 @@ bruteForceBracket(const std::function& eq, eq_low = eq_high; } if (bracket_found) { - // std::cout << "BFB: low: " << low << " high: " << high << std::endl; - // std::cout << "BFB: low_value: " << eq_low << " high_value: " << eq_high << std::endl; - // double low_value = eq(low); - // double high_value = eq(high); - // std::cout << "CHK: low_value: " << low_value<< " high_value: " << high_value << std::endl; deferred_logger.debug( " brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) + " high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high)); @@ -1010,6 +1011,47 @@ getFloIPR(const WellState& well_state, detail::getFlo(table, aqua_b, liquid_b, vapour_b)); } +bool +WellBhpThpCalculator:: +bruteForceBracketCommonTHP(const std::function& eq, + const std::array& range, + double& low, double& high, + std::optional& approximate_solution, + const double& limit, + DeferredLogger& deferred_logger) +{ + bool bracket_found = false; + // bool approximate_solution_found = false; + low = range[0]; + high = range[1]; + const int sample_number = 300;//10000; //300; //10000; + const double interval = (high - low) / sample_number; + double eq_low = eq(low); + double eq_high = 0.0; + for (int i = 0; i < sample_number + 1; ++i) { + high = range[0] + interval * i; + eq_high = eq(high); + if ( (std::fabs(eq_high) < limit)) { + // approximate_solution_found = true; + approximate_solution = high; + break; + } + if (eq_high * eq_low <= 0.) { + bracket_found = true; + break; + } + low = high; + eq_low = eq_high; + } + + if (bracket_found) { + deferred_logger.debug( + " brute force solve found low " + std::to_string(low) + " with eq_low " + std::to_string(eq_low) + + " high " + std::to_string(high) + " with eq_high " + std::to_string(eq_high)); + } + return bracket_found; +} + template class WellBhpThpCalculator; #define INSTANCE(...) \ diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index 24806d31e..72ca430a1 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -117,18 +117,21 @@ public: std::pair getFloIPR(const WellState& well_state, const Well& well, - const SummaryState& summary_state) const; - //! \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); + const SummaryState& summary_state) const; - //! \brief Find limits using brute-force solver. - static bool bruteForceBracketHigh(const std::function& eq, - const std::array& range, - double& low, double& high, - DeferredLogger& deferred_logger); + //! \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); private: //! \brief Compute BHP from THP limit for an injector - implementation. diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 6ca7bd441..1c1cc40a2 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -935,7 +935,7 @@ computeNetworkPressures(const Network::ExtNetwork& network, OpmLog::debug(oss.str()); #endif } else { - // Table number specified as 9999 in the deck, no pressure loss. + // Table number specified as 9999 in the deck, no pressure loss. if (network.node(node).as_choke()){ // Node pressure is set to the common THP of the wells. // The choke pressure must be non-negative therefore the node pressure of diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index bd78bdcfd..1add5ce2e 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -382,7 +382,18 @@ public: const WellProductionControls& prod_controls, WellState& well_state, const GroupState& group_state, - DeferredLogger& deferred_logger) = 0; + DeferredLogger& deferred_logger, + const bool fixed_control = false, + const bool fixed_status = false) = 0; + + bool solveWellWithTHPConstraint(const Simulator& simulator, + const double dt, + const Well::InjectionControls& inj_controls, + const Well::ProductionControls& prod_controls, + WellState& well_state, + const GroupState& group_state, + DeferredLogger& deferred_logger); + protected: // simulation parameters const ModelParameters& param_; diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 3e6464e11..6ab94cfda 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -487,6 +487,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/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 07049d76b..ce484b96b 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -103,6 +103,7 @@ public: void setWsolvent(const Scalar wsolvent); void setDynamicThpLimit(const Scalar thp_limit); std::optional getDynamicThpLimit() const; + void setDynamicThpLimit(const std::optional thp_limit); void updatePerforatedCell(std::vector& is_cell_perforated); /// Returns true if the well has one or more THP limits/constraints. From 393c70a83e04d7d8864d0356e445abf79ea5d755 Mon Sep 17 00:00:00 2001 From: Paul Date: Mon, 19 Feb 2024 22:09:33 +0100 Subject: [PATCH 6/8] clean up and improvements according reviewer comments --- opm/simulators/wells/BlackoilWellModel.hpp | 2 +- .../wells/BlackoilWellModelGeneric.cpp | 1 + .../wells/BlackoilWellModel_impl.hpp | 111 +++++------------- opm/simulators/wells/WellBhpThpCalculator.cpp | 4 +- opm/simulators/wells/WellGroupHelpers.cpp | 2 +- opm/simulators/wells/WellInterface_impl.hpp | 4 +- 6 files changed, 35 insertions(+), 89 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index b65004ea8..703c5aad5 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -413,7 +413,7 @@ template class WellContributions; Scalar gravity_{}; std::vector depth_{}; bool alternative_well_rate_init_{}; - double well_group_thp_prev_{0.0}; + std::optional well_group_thp_calc_; std::unique_ptr rateConverter_{}; std::map> regionalAveragePressureCalculator_{}; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 6f6eb88d3..32b9ea5ef 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1363,6 +1363,7 @@ updateNetworkPressures(const int reportStepIdx) } 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 2468dc96a..aa4714966 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1181,7 +1181,6 @@ namespace Opm { if (this->terminal_output_ && (network_update_iteration == iteration_to_relax) ) { local_deferredLogger.info(" we begin using relaxed tolerance for network update now after " + std::to_string(iteration_to_relax) + " iterations "); } - const bool relax_network_balance = network_update_iteration >= iteration_to_relax; std::tie(do_network_update, well_group_control_changed) = updateWellControlsAndNetworkIteration(mandatory_network_balance, relax_network_balance, dt,local_deferredLogger); @@ -1274,18 +1273,19 @@ namespace Opm { if (has_choke) { const auto& summary_state = this->ebosSimulator_.vanguard().summaryState(); const Group& group = schedule().getGroup(nodeName, reportStepIdx); - std::optional ctrl = group.productionControls(summary_state); - const auto cmode = ctrl->cmode; + const auto ctrl = group.productionControls(summary_state); + const auto cmode = ctrl.cmode; const auto pu = this->phase_usage_; - //PJPE: conversion factor res<-> surface rates TODO: to be calculated + //TODO: Auto choke combined with RESV control is not supported std::vector resv_coeff(pu.num_phases, 1.0); - // rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS. - double gratTargetFromSales = 0.0; //PJPE: check this + double gratTargetFromSales = 0.0; + if (group_state.has_grat_sales_target(group.name())) + gratTargetFromSales = group_state.grat_sales_target(group.name()); + WellGroupHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff, gratTargetFromSales, nodeName, group_state, group.has_gpmaint_control(cmode)); - double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); - //PJPE: TODO: include efficiency factor + const double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); auto mismatch = [&] (auto well_group_thp) { double well_group_rate(0.0); @@ -1294,32 +1294,18 @@ namespace Opm { std::string well_name = well->name(); auto& ws = well_state.well(well_name); if (group.hasWell(well_name)) { - rate = -tcalc.calcModeRateFromRates(ws.surface_rates); - std::cout << "(1) well: " << well_name << ", rate: " << rate << ", thp: " << ws.thp << ", Cmode: " << ws.production_cmode << std::endl; well->setDynamicThpLimit(well_group_thp); - const Well& well_ecl = wells_ecl_[well->indexOfWell()]; const auto inj_controls = Well::InjectionControls(0); const auto prod_controls = well_ecl.productionControls(summary_state); well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false); - // May be used instead if the wells of the sub-sea manifold has only THP constraints - // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); - rate = -tcalc.calcModeRateFromRates(ws.surface_rates); - if (ws.production_cmode == Opm::WellProducerCMode::THP) - ws.thp = well_group_thp; - - std::cout << "(2) well: " << well_name << ", rate: " << rate << ", thp: " << ws.thp << ", Cmode: " << ws.production_cmode << std::endl; well_group_rate += rate; } } - double diff = well_group_rate - orig_target; // (well_group_rate - orig_target)/orig_target; - std::cout << "well_group_thp: " << well_group_thp << " mismatch: " << diff << "\n" << std::endl; - return diff; + return well_group_rate - orig_target; }; - double well_group_rate(0.0); - double rate(0.0); double thp(0.0); double min_thp(1.0E8); double max_thp(0.0); @@ -1327,8 +1313,6 @@ namespace Opm { std::string well_name = well->name(); if (group.hasWell(well_name)) { auto& ws = well_state.well(well_name); - rate = -tcalc.calcModeRateFromRates(ws.surface_rates); - well_group_rate += rate; thp = ws.thp; min_thp = std::min(min_thp, thp); max_thp = std::max(max_thp, thp); @@ -1340,83 +1324,46 @@ namespace Opm { const double nodal_pressure = it->second; double well_group_thp = nodal_pressure; - if ((reportStepIdx <= 1) || (well_group_rate > 0.9*orig_target)) { - std::array range; // = {1.0E5, 150.0E5}; //PJPE what lower/upper bound to be taken here? - if (well_group_thp_prev_ == 0) { - range = {0.9*min_thp, 1.1*max_thp}; - } - else { - range = {0.9* this->well_group_thp_prev_, 1.1*this->well_group_thp_prev_}; - } + if (!this->well_group_thp_calc_.has_value() || this->well_group_thp_calc_ > nodal_pressure) { + std::array range; + this->well_group_thp_calc_.has_value() ? + range = {0.9 * this->well_group_thp_calc_.value(), 1.1 * this->well_group_thp_calc_.value()} : range = {0.9 * min_thp, 1.1 * max_thp}; double low, high; - // For testing - // low = 3.42767e+06;// 3.32833e+06; - // // high = 3.42767e+06;// 3.32833e+06; - // // std::cout << "low: " << low << " high: " << high << std::endl; - // double low_value_test = mismatch(low); - // low_value_test = mismatch(low); - // low_value_test = mismatch(low); - // std::cout << "low_value: " << low_value_test << " high_value: " << high_value_test << std::endl; - std::optional approximate_solution; - const double limit = 1.0E-3;// 1.0E-4; - local_deferredLogger.debug(" Using brute force search to bracket the common THP "); - std::cout << " Using brute force search to bracket the common THP " << std::endl; - const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, limit, local_deferredLogger); - - double low_value = mismatch(low); - double high_value = mismatch(high); + const double tolerance1 = 1.0E-3; + local_deferredLogger.debug("Using brute force search to bracket the common THP"); + const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger); if (approximate_solution.has_value()) { - well_group_thp = *approximate_solution; - double check = mismatch(well_group_thp); - local_deferredLogger.debug("Approximate value found: mismatch = " + std::to_string(check) + ", well_group_thp = " + std::to_string(well_group_thp)); + this->well_group_thp_calc_ = *approximate_solution; + local_deferredLogger.debug("Approximate common THP value found: " + std::to_string(this->well_group_thp_calc_.value())); } else if (finding_bracket) { - // PJPE: TODO what settings to take here? - std::cout << "Found bracket: [" << low << ", " << high << "]; mismatch: [" << low_value << ", "<< high_value << "]" << std::endl; - const double tolerance = 1.0E-3; + const double tolerance2 = 1.0E-3; const int max_iteration_solve = 100; int iteration = 0; - well_group_thp = RegulaFalsiBisection:: - solve(mismatch, low, high, max_iteration_solve, tolerance, iteration); - local_deferredLogger.debug( " bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " + - "mismatch = [" + std::to_string(low_value) + ", " + std::to_string(high_value) + "], " + - "iteration = " + std::to_string(iteration)); - double check = mismatch(well_group_thp); - local_deferredLogger.debug(" mismatch = " + std::to_string(check) + ", well_group_thp = " + std::to_string(well_group_thp)); + this->well_group_thp_calc_= RegulaFalsiBisection:: + solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration); + local_deferredLogger.debug(" bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " + + "iteration = " + std::to_string(iteration)); + local_deferredLogger.debug("Common THP value = " + std::to_string(this->well_group_thp_calc_.value())); } else { local_deferredLogger.warning("Common THP solve failed due to bracketing failure"); - well_group_thp = nodal_pressure; } - - std::cout << "well_group_thp = " << well_group_thp << ", nodal_pressure = " << nodal_pressure << std::endl; - - } else { - well_group_thp = nodal_pressure; } - well_group_thp = std::max(well_group_thp, nodal_pressure); + well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure); for (auto& well : this->well_container_) { std::string well_name = well->name(); if (group.hasWell(well_name)) { auto& ws = well_state.well(well_name); - std::cout << "Cmode of Well: " << well_name << " is " << ws.production_cmode <setDynamicThpLimit(well_group_thp); ws.thp = well_group_thp; } - // const Well& well_ecl = wells_ecl_[well->indexOfWell()]; - // const auto inj_controls = Well::InjectionControls(0); - // const auto prod_controls = well_ecl.productionControls(summary_state); - // // well->updateWellStateWithTHPTargetProd(this->ebosSimulator_, well_state, local_deferredLogger); - // well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state_copy, group_state, local_deferredLogger); - if (ws.production_cmode == Opm::WellProducerCMode::THP) - ws.thp = well_group_thp; } } - well_group_thp_prev_ = well_group_thp; // Use the common THP in computeNetworkPressures group_state.update_well_group_thp(nodeName, well_group_thp); @@ -2555,10 +2502,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/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 9017bedb2..1c501ee18 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -1021,10 +1021,9 @@ bruteForceBracketCommonTHP(const std::function& eq, DeferredLogger& deferred_logger) { bool bracket_found = false; - // bool approximate_solution_found = false; low = range[0]; high = range[1]; - const int sample_number = 300;//10000; //300; //10000; + const int sample_number = 300; const double interval = (high - low) / sample_number; double eq_low = eq(low); double eq_high = 0.0; @@ -1032,7 +1031,6 @@ bruteForceBracketCommonTHP(const std::function& eq, high = range[0] + interval * i; eq_high = eq(high); if ( (std::fabs(eq_high) < limit)) { - // approximate_solution_found = true; approximate_solution = high; break; } diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 1c1cc40a2..6ca7bd441 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -935,7 +935,7 @@ computeNetworkPressures(const Network::ExtNetwork& network, OpmLog::debug(oss.str()); #endif } else { - // Table number specified as 9999 in the deck, no pressure loss. + // Table number specified as 9999 in the deck, no pressure loss. if (network.node(node).as_choke()){ // Node pressure is set to the common THP of the wells. // The choke pressure must be non-negative therefore the node pressure of diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 098c6ff1c..ec451bc56 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1416,10 +1416,10 @@ 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 6e76602e8f1fa1f0eeea34dd57379539edf6eca0 Mon Sep 17 00:00:00 2001 From: Paul Date: Sun, 24 Mar 2024 20:47:46 +0100 Subject: [PATCH 7/8] changed assessing safe THP range --- opm/simulators/wells/BlackoilWellModel.hpp | 2 +- .../wells/BlackoilWellModel_impl.hpp | 94 +++++++++++-------- opm/simulators/wells/GroupState.cpp | 5 +- opm/simulators/wells/WellBhpThpCalculator.cpp | 49 +++++++--- opm/simulators/wells/WellBhpThpCalculator.hpp | 26 +++-- opm/simulators/wells/WellInterface.hpp | 25 +---- 6 files changed, 106 insertions(+), 95 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 703c5aad5..92afed768 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -466,7 +466,7 @@ template class WellContributions; const double dt, DeferredLogger& local_deferredLogger); - void computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger); + void 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 aa4714966..23c405a9b 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -49,12 +49,9 @@ #include #include -<<<<<<< HEAD #if COMPILE_BDA_BRIDGE #include #endif -======= ->>>>>>> f92db77d3 (moved common thp calculation to updateWellControls) #if HAVE_MPI #include @@ -1258,8 +1255,10 @@ namespace Opm { BlackoilWellModel:: computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger) { - const int reportStepIdx = this->ebosSimulator_.episodeIndex(); - const auto& network = schedule()[reportStepIdx].network(); + const int reportStepIdx = this->simulator_.episodeIndex(); + const auto& network = this->schedule()[reportStepIdx].network(); + const auto& balance = this->schedule()[reportStepIdx].network_balance(); + const Scalar thp_tolerance = balance.thp_tolerance(); if (!network.active()) { return; @@ -1271,75 +1270,87 @@ namespace Opm { for (const std::string& nodeName : network.node_names()) { const bool has_choke = network.node(nodeName).as_choke(); if (has_choke) { - const auto& summary_state = this->ebosSimulator_.vanguard().summaryState(); - const Group& group = schedule().getGroup(nodeName, reportStepIdx); + 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_; //TODO: Auto choke combined with RESV control is not supported - std::vector resv_coeff(pu.num_phases, 1.0); - double gratTargetFromSales = 0.0; + 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()); - WellGroupHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff, + WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff, gratTargetFromSales, nodeName, group_state, group.has_gpmaint_control(cmode)); - const double orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); + const Scalar orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); auto mismatch = [&] (auto well_group_thp) { - double well_group_rate(0.0); - double rate(0.0); + Scalar well_group_rate(0.0); + Scalar rate(0.0); for (auto& well : this->well_container_) { std::string well_name = well->name(); auto& ws = well_state.well(well_name); if (group.hasWell(well_name)) { well->setDynamicThpLimit(well_group_thp); - const Well& well_ecl = wells_ecl_[well->indexOfWell()]; + const Well& well_ecl = this->wells_ecl_[well->indexOfWell()]; const auto inj_controls = Well::InjectionControls(0); const auto prod_controls = well_ecl.productionControls(summary_state); - well->iterateWellEqWithSwitching(this->ebosSimulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false); + well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false); rate = -tcalc.calcModeRateFromRates(ws.surface_rates); well_group_rate += rate; } } - return well_group_rate - orig_target; + return (well_group_rate - orig_target)/orig_target; }; - double thp(0.0); - double min_thp(1.0E8); - double max_thp(0.0); - for (auto& well : this->well_container_) { - std::string well_name = well->name(); - if (group.hasWell(well_name)) { - auto& ws = well_state.well(well_name); - thp = ws.thp; - min_thp = std::min(min_thp, thp); - max_thp = std::max(max_thp, thp); + Scalar 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(); + WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp); + // Narrow down the bracket + Scalar 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 = node_pressures_.find((*upbranch).uptree_node()); - const double nodal_pressure = it->second; - double well_group_thp = nodal_pressure; + const auto it = this->node_pressures_.find((*upbranch).uptree_node()); + const Scalar nodal_pressure = it->second; + Scalar well_group_thp = nodal_pressure; if (!this->well_group_thp_calc_.has_value() || this->well_group_thp_calc_ > nodal_pressure) { - std::array range; + // The bracket is based on the initial bracket or on a range based on a previous calculated common group thp + std::array range; this->well_group_thp_calc_.has_value() ? - range = {0.9 * this->well_group_thp_calc_.value(), 1.1 * this->well_group_thp_calc_.value()} : range = {0.9 * min_thp, 1.1 * max_thp}; + range = {0.9 * this->well_group_thp_calc_.value(), 1.1 * this->well_group_thp_calc_.value()} : + range = range_initial; - double low, high; - std::optional approximate_solution; - const double tolerance1 = 1.0E-3; + Scalar low, high; + std::optional approximate_solution; + const Scalar tolerance1 = thp_tolerance; local_deferredLogger.debug("Using brute force search to bracket the common THP"); - const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger); + const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger); if (approximate_solution.has_value()) { this->well_group_thp_calc_ = *approximate_solution; local_deferredLogger.debug("Approximate common THP value found: " + std::to_string(this->well_group_thp_calc_.value())); } else if (finding_bracket) { - const double tolerance2 = 1.0E-3; + const Scalar tolerance2 = thp_tolerance; const int max_iteration_solve = 100; int iteration = 0; this->well_group_thp_calc_= RegulaFalsiBisection:: @@ -1348,11 +1359,12 @@ namespace Opm { "iteration = " + std::to_string(iteration)); local_deferredLogger.debug("Common THP value = " + std::to_string(this->well_group_thp_calc_.value())); } else { - local_deferredLogger.warning("Common THP solve failed due to bracketing failure"); + this->well_group_thp_calc_ = {}; + local_deferredLogger.debug("Common THP solve failed due to bracketing failure"); } } - - well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure); + this->well_group_thp_calc_.has_value() ? + well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure) : well_group_thp = nodal_pressure; for (auto& well : this->well_container_) { std::string well_name = well->name(); @@ -1365,7 +1377,7 @@ namespace Opm { } } - // Use the common THP in computeNetworkPressures + // Use the common group THP in computeNetworkPressures group_state.update_well_group_thp(nodeName, well_group_thp); } } @@ -2031,7 +2043,7 @@ namespace Opm { // network related bool more_network_update = false; if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) { - const double dt = this->ebosSimulator_.timeStepSize(); + 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); const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx); diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index a46ed4ef1..ee3fbf5f3 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -101,8 +101,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; @@ -119,6 +120,8 @@ GroupState::well_group_thp(const std::string& gname) const return group_iter->second; } +//------------------------------------------------------------------------- + template bool GroupState:: has_production_reduction_rates(const std::string& gname) const diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 1c501ee18..24377c096 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -875,15 +875,9 @@ bruteForceBracket(const std::function& eq, low = range[0]; high = range[1]; const int sample_number = 200; -<<<<<<< HEAD const Scalar interval = (high - low) / sample_number; Scalar eq_low = eq(low); Scalar eq_high = 0.0; -======= - const double interval = (high - low) / sample_number; - double eq_low = eq(low); - double eq_high = 0.0; ->>>>>>> 3ce8d7eec (moved common thp calculation to updateWellControls) for (int i = 0; i < sample_number + 1; ++i) { high = range[0] + interval * i; eq_high = eq(high); @@ -1011,22 +1005,23 @@ getFloIPR(const WellState& well_state, detail::getFlo(table, aqua_b, liquid_b, vapour_b)); } +template bool -WellBhpThpCalculator:: -bruteForceBracketCommonTHP(const std::function& eq, - const std::array& range, - double& low, double& high, - std::optional& approximate_solution, - const double& limit, +WellBhpThpCalculator:: +bruteForceBracketCommonTHP(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& high, + std::optional& approximate_solution, + const Scalar& limit, DeferredLogger& deferred_logger) { bool bracket_found = false; low = range[0]; high = range[1]; const int sample_number = 300; - const double interval = (high - low) / sample_number; - double eq_low = eq(low); - double eq_high = 0.0; + const Scalar interval = (high - low) / sample_number; + Scalar eq_low = eq(low); + Scalar eq_high = 0.0; for (int i = 0; i < sample_number + 1; ++i) { high = range[0] + interval * i; eq_high = eq(high); @@ -1050,6 +1045,30 @@ bruteForceBracketCommonTHP(const std::function& eq, return bracket_found; } +template +bool +WellBhpThpCalculator:: +bruteForceBracketCommonTHP(const std::function& eq, + Scalar& min_thp, Scalar& max_thp) +{ + bool bracket_found = false; + const int sample_number = 1000; + const 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); + if (eq_high * eq_low <= 0.) { + bracket_found = true; + min_thp = max_thp - interval; + break; + } + eq_low = eq_high; + } + return bracket_found; +} + template class WellBhpThpCalculator; #define INSTANCE(...) \ diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index 72ca430a1..057d19ff0 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -120,19 +120,23 @@ public: const SummaryState& summary_state) const; //! \brief Find limits using brute-force solver. - static bool bruteForceBracket(const std::function& eq, - const std::array& range, - double& low, double& high, + static bool bruteForceBracket(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& 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, + static bool bruteForceBracketCommonTHP(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& high, + std::optional& approximate_solution, + const Scalar& limit, DeferredLogger& deferred_logger); + //! \brief Find limits using brute-force solver. + static bool bruteForceBracketCommonTHP(const std::function& eq, + Scalar& min_thp, Scalar& max_thp); + private: //! \brief Compute BHP from THP limit for an injector - implementation. template @@ -169,12 +173,6 @@ private: std::optional& approximate_solution, DeferredLogger& deferred_logger) const; - //! \brief Find limits using brute-force solver. - static bool bruteForceBracket(const std::function& eq, - const std::array& range, - Scalar& low, Scalar& high, - DeferredLogger& deferred_logger); - Scalar findThpFromBhpIteratively(const std::function& thp_func, const Scalar bhp, const Scalar thp_limit, diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 1add5ce2e..66e6950f7 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -374,26 +374,15 @@ public: void updateConnectionTransmissibilityFactor(const Simulator& simulator, SingleWellState& ws) const; - SingleWellState& ws) const; - virtual bool iterateWellEqWithSwitching(const Simulator& simulator, const double dt, const WellInjectionControls& inj_controls, const WellProductionControls& prod_controls, - WellState& well_state, - const GroupState& group_state, + WellState& well_state, + const GroupState& group_state, DeferredLogger& deferred_logger, const bool fixed_control = false, const bool fixed_status = false) = 0; - - bool solveWellWithTHPConstraint(const Simulator& simulator, - const double dt, - const Well::InjectionControls& inj_controls, - const Well::ProductionControls& prod_controls, - WellState& well_state, - const GroupState& group_state, - DeferredLogger& deferred_logger); - protected: // simulation parameters const ModelParameters& param_; @@ -448,16 +437,6 @@ protected: const GroupState& group_state, DeferredLogger& deferred_logger) = 0; - virtual bool iterateWellEqWithSwitching(const Simulator& simulator, - const double dt, - const WellInjectionControls& inj_controls, - const WellProductionControls& prod_controls, - WellState& well_state, - const GroupState& group_state, - DeferredLogger& deferred_logger, - const bool fixed_control = false, - const bool fixed_status = false) = 0; - virtual void updateIPRImplicit(const Simulator& simulator, WellState& well_state, DeferredLogger& deferred_logger) = 0; From f9d82c60427e3fdc5b0cdf359376e1df2f622e03 Mon Sep 17 00:00:00 2001 From: Paul Date: Tue, 25 Jun 2024 11:41:34 +0200 Subject: [PATCH 8/8] adressing several reviewers comments --- opm/simulators/wells/BlackoilWellModel.hpp | 3 +- .../wells/BlackoilWellModel_impl.hpp | 77 ++++++++++--------- opm/simulators/wells/GroupState.cpp | 2 +- opm/simulators/wells/GroupState.hpp | 4 +- opm/simulators/wells/WellBhpThpCalculator.cpp | 4 +- opm/simulators/wells/WellBhpThpCalculator.hpp | 13 ++-- opm/simulators/wells/WellGroupHelpers.cpp | 7 +- opm/simulators/wells/WellInterfaceGeneric.cpp | 7 -- 8 files changed, 54 insertions(+), 63 deletions(-) diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 92afed768..f877a8c80 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -413,8 +413,7 @@ template class WellContributions; Scalar gravity_{}; std::vector depth_{}; bool alternative_well_rate_init_{}; - std::optional well_group_thp_calc_; - + std::map well_group_thp_calc_; std::unique_ptr rateConverter_{}; std::map> regionalAveragePressureCalculator_{}; diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 23c405a9b..d00b31b3b 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -1282,40 +1282,49 @@ namespace Opm { gratTargetFromSales = group_state.grat_sales_target(group.name()); WGHelpers::TargetCalculator tcalc(cmode, pu, resv_coeff, - gratTargetFromSales, nodeName, group_state, - group.has_gpmaint_control(cmode)); + gratTargetFromSales, nodeName, group_state, + group.has_gpmaint_control(cmode)); const Scalar orig_target = tcalc.groupTarget(ctrl, local_deferredLogger); - auto mismatch = [&] (auto well_group_thp) { - Scalar well_group_rate(0.0); + auto mismatch = [&] (auto group_thp) { + Scalar group_rate(0.0); Scalar rate(0.0); for (auto& well : this->well_container_) { std::string well_name = well->name(); auto& ws = well_state.well(well_name); if (group.hasWell(well_name)) { - well->setDynamicThpLimit(well_group_thp); + well->setDynamicThpLimit(group_thp); const Well& well_ecl = this->wells_ecl_[well->indexOfWell()]; const auto inj_controls = Well::InjectionControls(0); const auto prod_controls = well_ecl.productionControls(summary_state); well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false); rate = -tcalc.calcModeRateFromRates(ws.surface_rates); - well_group_rate += rate; + group_rate += rate; } } - return (well_group_rate - orig_target)/orig_target; + return (group_rate - orig_target)/orig_target; }; - Scalar min_thp, max_thp; - std::array range_initial; + const auto upbranch = network.uptree_branch(nodeName); + const auto it = this->node_pressures_.find((*upbranch).uptree_node()); + const Scalar nodal_pressure = it->second; + Scalar well_group_thp = nodal_pressure; + + std::optional autochoke_thp; + if (auto iter = this->well_group_thp_calc_.find(nodeName); iter != this->well_group_thp_calc_.end()) { + autochoke_thp = this->well_group_thp_calc_.at(nodeName); + } + //Find an initial bracket - if (!this->well_group_thp_calc_.has_value()){ + std::array range_initial; + if (!autochoke_thp.has_value()){ + Scalar min_thp, max_thp; // 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(); WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp); // Narrow down the bracket @@ -1328,56 +1337,48 @@ namespace Opm { 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; - Scalar well_group_thp = nodal_pressure; - - if (!this->well_group_thp_calc_.has_value() || this->well_group_thp_calc_ > nodal_pressure) { - // The bracket is based on the initial bracket or on a range based on a previous calculated common group thp - std::array range; - this->well_group_thp_calc_.has_value() ? - range = {0.9 * this->well_group_thp_calc_.value(), 1.1 * this->well_group_thp_calc_.value()} : - range = range_initial; - + if (!autochoke_thp.has_value() || autochoke_thp.value() > nodal_pressure) { + // The bracket is based on the initial bracket or on a range based on a previous calculated group thp + std::array range = autochoke_thp.has_value() ? + std::array{0.9 * autochoke_thp.value(), 1.1 * autochoke_thp.value()} : range_initial; Scalar low, high; std::optional approximate_solution; const Scalar tolerance1 = thp_tolerance; - local_deferredLogger.debug("Using brute force search to bracket the common THP"); + local_deferredLogger.debug("Using brute force search to bracket the group THP"); const bool finding_bracket = WellBhpThpCalculator::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger); if (approximate_solution.has_value()) { - this->well_group_thp_calc_ = *approximate_solution; - local_deferredLogger.debug("Approximate common THP value found: " + std::to_string(this->well_group_thp_calc_.value())); + autochoke_thp = *approximate_solution; + local_deferredLogger.debug("Approximate group THP value found: " + std::to_string(autochoke_thp.value())); } else if (finding_bracket) { const Scalar tolerance2 = thp_tolerance; const int max_iteration_solve = 100; int iteration = 0; - this->well_group_thp_calc_= RegulaFalsiBisection:: + autochoke_thp = RegulaFalsiBisection:: solve(mismatch, low, high, max_iteration_solve, tolerance2, iteration); local_deferredLogger.debug(" bracket = [" + std::to_string(low) + ", " + std::to_string(high) + "], " + "iteration = " + std::to_string(iteration)); - local_deferredLogger.debug("Common THP value = " + std::to_string(this->well_group_thp_calc_.value())); + local_deferredLogger.debug("Group THP value = " + std::to_string(autochoke_thp.value())); } else { - this->well_group_thp_calc_ = {}; - local_deferredLogger.debug("Common THP solve failed due to bracketing failure"); + autochoke_thp.reset(); + local_deferredLogger.debug("Group THP solve failed due to bracketing failure"); } } - this->well_group_thp_calc_.has_value() ? - well_group_thp = std::max(this->well_group_thp_calc_.value(), nodal_pressure) : well_group_thp = nodal_pressure; + if (autochoke_thp.has_value()) { + well_group_thp_calc_[nodeName] = autochoke_thp.value(); + // Note: The node pressure of the auto-choke node is set to well_group_thp in computeNetworkPressures() + // and must be larger or equal to the pressure of the uptree node of its branch. + well_group_thp = std::max(autochoke_thp.value(), nodal_pressure); + } for (auto& well : this->well_container_) { std::string well_name = well->name(); if (group.hasWell(well_name)) { - auto& ws = well_state.well(well_name); - if (ws.production_cmode == Opm::WellProducerCMode::THP) { - well->setDynamicThpLimit(well_group_thp); - ws.thp = well_group_thp; - } + well->setDynamicThpLimit(well_group_thp); } } - // Use the common group THP in computeNetworkPressures + // Use the group THP in computeNetworkPressures(). group_state.update_well_group_thp(nodeName, well_group_thp); } } diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index ee3fbf5f3..a3e5c1101 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -110,7 +110,7 @@ GroupState::update_well_group_thp(const std::string& gname, const double& thp) } template -double GroupState:: +Scalar GroupState:: GroupState::well_group_thp(const std::string& gname) const { auto group_iter = this->group_thp.find(gname); diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index 97ab1d0f0..cc22cde8b 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -50,7 +50,7 @@ public: const std::vector& production_rates(const std::string& gname) const; void update_well_group_thp(const std::string& gname, const double& thp); - double well_group_thp(const std::string& gname) const; + Scalar well_group_thp(const std::string& gname) const; bool has_production_reduction_rates(const std::string& gname) const; void update_production_reduction_rates(const std::string& gname, @@ -203,7 +203,7 @@ private: std::map inj_vrep_rate; std::map m_grat_sales_target; std::map m_gpmaint_target; - std::map group_thp; + std::map group_thp; std::map, Group::InjectionCMode> injection_controls; WellContainer gpmaint_state; diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 24377c096..0ea2c1828 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -1052,8 +1052,8 @@ bruteForceBracketCommonTHP(const std::function& eq, Scalar& min_thp, Scalar& max_thp) { bool bracket_found = false; - const int sample_number = 1000; - const Scalar interval = 1E5; + 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) { diff --git a/opm/simulators/wells/WellBhpThpCalculator.hpp b/opm/simulators/wells/WellBhpThpCalculator.hpp index 057d19ff0..c1461260e 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.hpp +++ b/opm/simulators/wells/WellBhpThpCalculator.hpp @@ -119,12 +119,6 @@ public: const Well& well, const SummaryState& summary_state) const; - //! \brief Find limits using brute-force solver. - static bool bruteForceBracket(const std::function& eq, - const std::array& range, - Scalar& low, Scalar& high, - DeferredLogger& deferred_logger); - //! \brief Find limits using brute-force solver. static bool bruteForceBracketCommonTHP(const std::function& eq, const std::array& range, @@ -173,6 +167,13 @@ private: std::optional& approximate_solution, DeferredLogger& deferred_logger) const; + //! \brief Find limits using brute-force solver. + static bool bruteForceBracket(const std::function& eq, + const std::array& range, + Scalar& low, Scalar& high, + DeferredLogger& deferred_logger); + + Scalar findThpFromBhpIteratively(const std::function& thp_func, const Scalar bhp, const Scalar thp_limit, diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index 6ca7bd441..1ceb5a49a 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -937,11 +937,8 @@ computeNetworkPressures(const Network::ExtNetwork& network, } else { // Table number specified as 9999 in the deck, no pressure loss. if (network.node(node).as_choke()){ - // Node pressure is set to the common THP of the wells. - // The choke pressure must be non-negative therefore the node pressure of - // the auto-choke node must be larger or equal to the pressure of the uptree node of its branch - const auto group_thp = group_state.well_group_thp(node); - node_pressures[node] = group_thp >= up_press ? group_thp : up_press; + // Node pressure is set to the group THP. + node_pressures[node] = group_state.well_group_thp(node); } else { node_pressures[node] = up_press; } diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 6ab94cfda..3e6464e11 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -487,13 +487,6 @@ 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)