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