diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f8a47b0fad5..51fb806e0c1 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -208,6 +208,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/BlackoilWellModelGasLift.cpp opm/simulators/wells/BlackoilWellModelGeneric.cpp opm/simulators/wells/BlackoilWellModelGuideRates.cpp + opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp opm/simulators/wells/BlackoilWellModelNldd.cpp opm/simulators/wells/BlackoilWellModelRestart.cpp opm/simulators/wells/BlackoilWellModelWBP.cpp @@ -1131,6 +1132,9 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/BlackoilWellModelGasLift_impl.hpp opm/simulators/wells/BlackoilWellModelGeneric.hpp opm/simulators/wells/BlackoilWellModelGuideRates.hpp + opm/simulators/wells/BlackoilWellModelNetwork.hpp + opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp + opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp opm/simulators/wells/BlackoilWellModelNldd.hpp opm/simulators/wells/BlackoilWellModelNldd_impl.hpp opm/simulators/wells/BlackoilWellModelRestart.hpp diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index caf0fe2f660..04292545fe3 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -57,6 +57,7 @@ #include #include #include +#include #include #include #include @@ -281,10 +282,6 @@ template class WellContributions; bool updateWellControls(DeferredLogger& deferred_logger); - std::tuple - updateNetworks(const bool mandatory_network_balance, DeferredLogger& deferred_logger, const bool relax_network_tolerance = false); - - void updateAndCommunicate(const int reportStepIdx, const int iterationIdx, DeferredLogger& deferred_logger); @@ -425,6 +422,20 @@ template class WellContributions; setupRescoupScopedLogger(DeferredLogger& local_logger); #endif + + bool updateWellControlsAndNetwork(const bool mandatory_network_balance, + const double dt, + DeferredLogger& local_deferredLogger); + + // TODO: finding a better naming + void assembleWellEqWithoutIteration(const double dt, DeferredLogger& deferred_logger); + + const std::vector& B_avg() const + { return B_avg_; } + + const ModelParameters& param() const + { return param_; } + protected: Simulator& simulator_; @@ -458,7 +469,6 @@ template class WellContributions; Scalar gravity_{}; std::vector depth_{}; bool alternative_well_rate_init_{}; - std::map well_group_thp_calc_; std::unique_ptr rateConverter_{}; std::map> regionalAveragePressureCalculator_{}; @@ -468,9 +478,6 @@ template class WellContributions; // A flag to tell the convergence report whether we need to take another newton step bool network_needs_more_balancing_force_another_newton_iteration_{false}; - // Pre-step network solve at static reservoir conditions (group and well states might be updated) - void doPreStepNetworkRebalance(DeferredLogger& deferred_logger); - std::vector B_avg_{}; const EquilGrid& equilGrid() const @@ -494,12 +501,6 @@ template class WellContributions; const double dt, DeferredLogger& local_deferredLogger); - bool updateWellControlsAndNetwork(const bool mandatory_network_balance, - const double dt, - DeferredLogger& local_deferredLogger); - - bool computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger); - /// Update rank's notion of intersecting wells and their /// associate solution variables. /// @@ -549,9 +550,6 @@ template class WellContributions; void prepareWellsBeforeAssembling(const double dt, DeferredLogger& deferred_logger); - // TODO: finding a better naming - void assembleWellEqWithoutIteration(const double dt, DeferredLogger& deferred_logger); - void extractLegacyCellPvtRegionIndex_(); void extractLegacyDepth_(); @@ -574,6 +572,7 @@ template class WellContributions; private: BlackoilWellModelGasLift gaslift_; + BlackoilWellModelNetwork network_; BlackoilWellModelNldd* nldd_ = nullptr; //!< NLDD well model adapter (not owned) // These members are used to avoid reallocation in specific functions diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index bf1d174c4df..ceba6a2f647 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -58,6 +58,7 @@ #include #include #include +#include #include #include #include @@ -92,6 +93,7 @@ template BlackoilWellModelGeneric:: BlackoilWellModelGeneric(Schedule& schedule, BlackoilWellModelGasLiftGeneric& gaslift, + BlackoilWellModelNetworkGeneric& network, const SummaryState& summaryState, const EclipseState& eclState, const PhaseUsageInfo& pu, @@ -109,7 +111,13 @@ BlackoilWellModelGeneric(Schedule& schedule, , active_wgstate_(pu) , last_valid_wgstate_(pu) , nupcol_wgstate_(pu) - , group_state_helper_(this->wellState(), this->groupState(), this->schedule(), summaryState, guideRate_, pu) + , group_state_helper_(this->wellState(), + this->groupState(), + this->schedule(), + summaryState, + guideRate_, + pu) + , genNetwork_(network) { const auto numProcs = comm_.size(); @@ -127,17 +135,6 @@ BlackoilWellModelGeneric(Schedule& schedule, return (candidate == this->parallel_well_info_.end()) || (*candidate != value); }; - - const auto& node_pressures = eclState.getRestartNetworkPressures(); - if (node_pressures.has_value()) { - if constexpr (std::is_same_v) { - this->node_pressures_ = node_pressures.value(); - } else { - for (const auto& it : node_pressures.value()) { - this->node_pressures_[it.first] = it.second; - } - } - } } template @@ -186,13 +183,6 @@ wellsActive() const return wells_active_; } -template -bool BlackoilWellModelGeneric:: -networkActive() const -{ - return network_active_; -} - template bool BlackoilWellModelGeneric:: anyMSWellOpenLocal() const @@ -1197,51 +1187,10 @@ assignGroupValues(const int reportStepIdx, template void BlackoilWellModelGeneric:: -assignNodeValues(std::map& nodevalues, - const int reportStepIdx) const +commitWGState() { - nodevalues.clear(); - if (reportStepIdx < 0) return; - - for (const auto& [node, pressure] : node_pressures_) { - nodevalues.emplace(node, data::NodeData{pressure}); - // Assign node values of well groups to GPR:WELLNAME - const auto& sched = schedule(); - if (!sched.hasGroup(node, reportStepIdx)) { - continue; - } - const auto& group = sched.getGroup(node, reportStepIdx); - for (const std::string& wellname : group.wells()) { - nodevalues.emplace(wellname, data::NodeData{pressure}); - } - } - - const auto& network = schedule()[reportStepIdx].network(); - if (!network.active()) { - return; - } - - auto converged_pressures = this->groupStateHelper().computeNetworkPressures( - network, - *(this->vfp_properties_->getProd()), - this->comm_ - ); - for (const auto& [node, converged_pressure] : converged_pressures) { - auto it = nodevalues.find(node); - assert(it != nodevalues.end() ); - it->second.converged_pressure = converged_pressure; - // Assign node values of group to GPR:WELLNAME - const auto& sched = schedule(); - if (!sched.hasGroup(node, reportStepIdx)) { - continue; - } - const auto& group = sched.getGroup(node, reportStepIdx); - for (const std::string& wellname : group.wells()) { - auto it2 = nodevalues.find(wellname); - assert(it2 != nodevalues.end()); - it2->second.converged_pressure = converged_pressure; - } - } + this->last_valid_wgstate_ = this->active_wgstate_; + this->genNetwork_.commitState(); } template @@ -1252,7 +1201,7 @@ groupAndNetworkData(const int reportStepIdx) const auto grp_nwrk_values = data::GroupAndNetworkValues{}; this->assignGroupValues(reportStepIdx, grp_nwrk_values.groupData); - this->assignNodeValues(grp_nwrk_values.nodeData, reportStepIdx - 1); // Schedule state info at previous step + this->genNetwork_.assignNodeValues(grp_nwrk_values.nodeData, reportStepIdx - 1); // Schedule state info at previous step return grp_nwrk_values; } @@ -1435,46 +1384,6 @@ updateAndCommunicateGroupData(const int reportStepIdx, } } -template -void BlackoilWellModelGeneric:: -updateNetworkActiveState(const int report_step) { - const auto& network = schedule()[report_step].network(); - if (!network.active()) { - this->network_active_ = false; - return; - } - - bool network_active = false; - for (const auto& well : well_container_generic_) { - const bool is_partof_network = network.has_node(well->wellEcl().groupName()); - const bool prediction_mode = well->wellEcl().predictionMode(); - if (is_partof_network && prediction_mode) { - network_active = true; - break; - } - } - this->network_active_ = comm_.max(network_active); -} - -template -bool BlackoilWellModelGeneric:: -needPreStepNetworkRebalance(const int report_step) const -{ - const auto& network = schedule()[report_step].network(); - bool network_rebalance_necessary = false; - for (const auto& well : well_container_generic_) { - const bool is_partof_network = network.has_node(well->wellEcl().groupName()); - // TODO: we might find more relevant events to be included here (including network change events?) - const auto& events = this->wellState().well(well->indexOfWell()).events; - if (is_partof_network && events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE)) { - network_rebalance_necessary = true; - break; - } - } - network_rebalance_necessary = comm_.max(network_rebalance_necessary); - return network_rebalance_necessary; -} - template bool BlackoilWellModelGeneric:: forceShutWellByName(const std::string& wellname, @@ -1550,80 +1459,6 @@ inferLocalShutWells() } } -template -Scalar -BlackoilWellModelGeneric:: -updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor, const Scalar upper_update_bound) -{ - OPM_TIMEFUNCTION(); - // Get the network and return if inactive (no wells in network at this time) - const auto& network = schedule()[reportStepIdx].network(); - if (!network.active()) { - return 0.0; - } - - const auto previous_node_pressures = node_pressures_; - - node_pressures_ = this->groupStateHelper().computeNetworkPressures( - network, *(vfp_properties_->getProd()), comm_); - - // here, the network imbalance is the difference between the previous nodal pressure and the new nodal pressure - Scalar network_imbalance = 0.; - if (!this->networkActive()) - return network_imbalance; - - if (!previous_node_pressures.empty()) { - for (const auto& [name, new_pressure]: node_pressures_) { - if (previous_node_pressures.count(name) <= 0) { - if (std::abs(new_pressure) > network_imbalance) { - network_imbalance = std::abs(new_pressure); - } - continue; - } - const auto pressure = previous_node_pressures.at(name); - const Scalar change = (new_pressure - pressure); - if (std::abs(change) > network_imbalance) { - network_imbalance = std::abs(change); - } - // We dampen the nodal pressure change during one iteration since our nodal pressure calculation - // is somewhat explicit. There is a relative dampening factor applied to the update value, and also - // the maximum update is limited (to 5 bar by default, can be changed with --network-max-pressure-update-in-bars). - const Scalar damped_change = std::min(damping_factor * std::abs(change), upper_update_bound); - const Scalar sign = change > 0 ? 1. : -1.; - node_pressures_[name] = pressure + sign * damped_change; - } - } else { - for (const auto& [name, pressure]: node_pressures_) { - if (std::abs(pressure) > network_imbalance) { - network_imbalance = std::abs(pressure); - } - } - } - - 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. - if (well->isProducer() && well->wellEcl().predictionMode()) { - const auto it = node_pressures_.find(well->wellEcl().groupName()); - if (it != node_pressures_.end()) { - // The well belongs to a group with has a network pressure constraint, - // set the dynamic THP constraint of the well accordingly. - const Scalar new_limit = it->second; - well->setDynamicThpLimit(new_limit); - SingleWellState& ws = this->wellState()[well->indexOfWell()]; - const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP; - // TODO: not sure why the thp is NOT updated properly elsewhere - if (thp_is_limit) { - ws.thp = well->getTHPConstraint(summaryState_); - } - } - } - } - return network_imbalance; -} - template void BlackoilWellModelGeneric:: calculateEfficiencyFactors(const int reportStepIdx) @@ -1801,31 +1636,6 @@ runWellPIScaling(const int reportStepIdx, this->last_run_wellpi_ = reportStepIdx; } -template -bool BlackoilWellModelGeneric:: -shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const -{ - // if network is not active, we do not need to balance the network - const auto& network = schedule()[reportStepIdx].network(); - if (!network.active()) { - return false; - } - - const auto& balance = schedule()[reportStepIdx].network_balance(); - if (balance.mode() == Network::Balance::CalcMode::TimeStepStart) { - return iterationIdx == 0; - } else if (balance.mode() == Network::Balance::CalcMode::NUPCOL) { - const int nupcol = schedule()[reportStepIdx].nupcol(); - return iterationIdx < nupcol; - } else { - // We do not support any other rebalancing modes, - // i.e. TimeInterval based rebalancing is not available. - // This should be warned about elsewhere, so we choose to - // avoid spamming with a warning here. - return false; - } -} - template std::vector BlackoilWellModelGeneric:: getCellsForConnections(const Well& well) const @@ -2123,8 +1933,7 @@ operator==(const BlackoilWellModelGeneric& rhs) const && this->last_run_wellpi_ == rhs.last_run_wellpi_ && this->local_shut_wells_ == rhs.local_shut_wells_ && this->closed_this_step_ == rhs.closed_this_step_ - && this->node_pressures_ == rhs.node_pressures_ - && this->last_valid_node_pressures_ == rhs.last_valid_node_pressures_ + && this->genNetwork_ == rhs.genNetwork_ && this->prev_inj_multipliers_ == rhs.prev_inj_multipliers_ && this->active_wgstate_ == rhs.active_wgstate_ && this->last_valid_wgstate_ == rhs.last_valid_wgstate_ diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index c9ecd581e3c..4260ae9c250 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -60,6 +60,7 @@ namespace Opm { class DeferredLogger; class EclipseState; template class BlackoilWellModelGasLiftGeneric; + template class BlackoilWellModelNetworkGeneric; template class GasLiftGroupInfo; template class GasLiftSingleWellGeneric; template class GasLiftWellState; @@ -97,6 +98,7 @@ class BlackoilWellModelGeneric public: BlackoilWellModelGeneric(Schedule& schedule, BlackoilWellModelGasLiftGeneric& gaslift, + BlackoilWellModelNetworkGeneric& network, const SummaryState& summaryState, const EclipseState& eclState, const PhaseUsageInfo& phase_usage, @@ -123,9 +125,6 @@ class BlackoilWellModelGeneric //! \brief Returns true if well is defined, open and has connections on current rank. bool hasOpenLocalWell(const std::string& well_name) const; - /// return true if network is active (at least one network well in prediction mode) - bool networkActive() const; - // whether there exists any multisegment well open on this process bool anyMSWellOpenLocal() const; @@ -208,22 +207,10 @@ class BlackoilWellModelGeneric with storeWellState() can then subsequently be recovered with the resetWellState() method. */ - void commitWGState() - { - this->last_valid_wgstate_ = this->active_wgstate_; - this->last_valid_node_pressures_ = this->node_pressures_; - } + void commitWGState(); data::GroupAndNetworkValues groupAndNetworkData(const int reportStepIdx) const; - /// Checks if network is active (at least one network well on prediction). - void updateNetworkActiveState(const int report_step); - - /// Checks if there are reasons to perform a pre-step network re-balance. - /// (Currently, the only reasons are network well status changes.) - /// (TODO: Consider if adding network change events would be helpful.) - bool needPreStepNetworkRebalance(const int report_step) const; - /// Shut down any single well /// Returns true if the well was actually found and shut. bool forceShutWellByName(const std::string& wellname, @@ -250,9 +237,6 @@ class BlackoilWellModelGeneric bool reportStepStarts() const { return report_step_starts_; } - bool shouldBalanceNetwork(const int reportStepIndex, - const int iterationIdx) const; - void updateClosedWellsThisStep(const std::string& well_name) const { this->closed_this_step_.insert(well_name); @@ -270,8 +254,7 @@ class BlackoilWellModelGeneric serializer(local_shut_wells_); serializer(closed_this_step_); serializer(guideRate_); - serializer(node_pressures_); - serializer(last_valid_node_pressures_); + serializer(genNetwork_); serializer(prev_inj_multipliers_); serializer(active_wgstate_); serializer(last_valid_wgstate_); @@ -306,6 +289,22 @@ class BlackoilWellModelGeneric GroupStateHelperType& groupStateHelper() { return group_state_helper_; } const GroupStateHelperType& groupStateHelper() const { return group_state_helper_; } + const VFPProperties& getVFPProperties() const + { + return *vfp_properties_; + } + + void updateAndCommunicateGroupData(const int reportStepIdx, + const int iterationIdx, + const Scalar tol_nupcol, + // we only want to update the wellgroup target + // after the groups have found their controls + const bool update_wellgrouptarget, + DeferredLogger& deferred_logger); + + const EclipseState& eclState() const + { return eclState_; } + protected: /* The dynamic state of the well model is maintained with an instance @@ -362,9 +361,9 @@ class BlackoilWellModelGeneric void resetWGState() { this->active_wgstate_ = this->last_valid_wgstate_; - this->node_pressures_ = this->last_valid_node_pressures_; // Update helper pointers to reference the restored active state this->group_state_helper_.updateState(this->wellState(), this->groupState()); + this->genNetwork_.resetState(); } /* @@ -389,10 +388,6 @@ class BlackoilWellModelGeneric bool wasDynamicallyShutThisTimeStep(const int well_index) const; - Scalar updateNetworkPressures(const int reportStepIdx, - const Scalar damping_factor, - const Scalar update_upper_bound); - void updateWsolvent(const Group& group, const int reportStepIdx, const WellState& wellState); @@ -435,8 +430,6 @@ class BlackoilWellModelGeneric data::GroupData& gdata) const; void assignGroupValues(const int reportStepIdx, std::map& gvalues) const; - void assignNodeValues(std::map& nodevalues, - const int reportStepIdx) const; void calculateEfficiencyFactors(const int reportStepIdx); @@ -456,12 +449,6 @@ class BlackoilWellModelGeneric const int max_number_of_group_switch, const bool update_group_switching_log); - void updateAndCommunicateGroupData(const int reportStepIdx, - const int iterationIdx, - const Scalar tol_nupcol, - const bool update_wellgrouptarget, // we only want to update the wellgrouptarget after the groups have found their controls - DeferredLogger& deferred_logger); - void inferLocalShutWells(); void setRepRadiusPerfLength(); @@ -529,11 +516,9 @@ class BlackoilWellModelGeneric BlackoilWellModelGasLiftGeneric& gen_gaslift_; BlackoilWellModelWBP wbp_; - const PhaseUsageInfo& phase_usage_info_; bool terminal_output_{false}; bool wells_active_{false}; - bool network_active_{false}; bool initial_step_{}; bool report_step_starts_{}; @@ -568,11 +553,6 @@ class BlackoilWellModelGeneric GuideRate guideRate_; std::unique_ptr> vfp_properties_{}; - // Network pressures for output and initialization - std::map node_pressures_; - // Valid network pressures for output and initialization for safe restart after failed iterations - std::map last_valid_node_pressures_; - // previous injection multiplier, it is used in the injection multiplier calculation for WINJMULT keyword std::unordered_map> prev_inj_multipliers_; @@ -599,6 +579,8 @@ class BlackoilWellModelGeneric // Store map of group name and close offending well for output std::map> closed_offending_wells_; + BlackoilWellModelNetworkGeneric& genNetwork_; + private: WellInterfaceGeneric* getGenWell(const std::string& well_name); diff --git a/opm/simulators/wells/BlackoilWellModelNetwork.hpp b/opm/simulators/wells/BlackoilWellModelNetwork.hpp new file mode 100644 index 00000000000..5984c769516 --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelNetwork.hpp @@ -0,0 +1,82 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BLACKOILWELLMODEL_NETWORK_HEADER_INCLUDED +#define OPM_BLACKOILWELLMODEL_NETWORK_HEADER_INCLUDED + +#include +#include + +#include + +#include +#include +#include + +namespace Opm { + class DeferredLogger; + template class BlackoilWellModel; +} + +namespace Opm { + +/// Class for handling the blackoil well network model. +template +class BlackoilWellModelNetwork : + public BlackoilWellModelNetworkGeneric, + typename GetPropType::IndexTraitsType> +{ + using BaseType = + BlackoilWellModelNetworkGeneric, + typename GetPropType::IndexTraitsType>; + + using FluidSystem = GetPropType; + using Indices = GetPropType; + using IndexTraits = typename FluidSystem::IndexTraitsType; + using Scalar = GetPropType; + +public: + BlackoilWellModelNetwork(BlackoilWellModel& well_model); + + std::tuple + update(const bool mandatory_network_balance, + DeferredLogger& deferred_logger, + const bool relax_network_tolerance = false); + + // Pre-step network solve at static reservoir conditions (group and well states might be updated) + void doPreStepRebalance(DeferredLogger& deferred_logger); + +protected: + /// 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 + bool computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger); + + BlackoilWellModel& well_model_; + std::map well_group_thp_calc_; +}; + +} // namespace Opm + +#include "BlackoilWellModelNetwork_impl.hpp" + +#endif // OPM_BLACKOILWELLMODEL_NETWORK_HEADER_INCLUDED diff --git a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp new file mode 100644 index 00000000000..bbaf9097aac --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.cpp @@ -0,0 +1,467 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include + +#include + +namespace Opm { + +template +BlackoilWellModelNetworkGeneric:: +BlackoilWellModelNetworkGeneric(BlackoilWellModelGeneric& well_model) + : well_model_(well_model) +{ + this->setFromRestart(well_model_.eclState().getRestartNetworkPressures()); +} + +template +void BlackoilWellModelNetworkGeneric:: +setFromRestart(const std::optional>& node_pressures) +{ + if (node_pressures.has_value()) { + if constexpr (std::is_same_v) { + this->node_pressures_ = node_pressures.value(); + } else { + for (const auto& it : node_pressures.value()) { + this->node_pressures_[it.first] = it.second; + } + } + } +} + +template +void BlackoilWellModelNetworkGeneric:: +updateActiveState(const int report_step) +{ + const auto& network = well_model_.schedule()[report_step].network(); + if (!network.active()) { + this->active_ = false; + return; + } + + bool network_active = false; + for (const auto& well : well_model_.genericWells()) { + const bool is_partof_network = network.has_node(well->wellEcl().groupName()); + const bool prediction_mode = well->wellEcl().predictionMode(); + if (is_partof_network && prediction_mode) { + network_active = true; + break; + } + } + this->active_ = well_model_.comm().max(network_active); +} + +template +bool BlackoilWellModelNetworkGeneric:: +needPreStepRebalance(const int report_step) const +{ + const auto& network = well_model_.schedule()[report_step].network(); + bool network_rebalance_necessary = false; + for (const auto& well : well_model_.genericWells()) { + const bool is_partof_network = network.has_node(well->wellEcl().groupName()); + // TODO: we might find more relevant events to be included here (including network change events?) + const auto& events = well_model_.wellState().well(well->indexOfWell()).events; + if (is_partof_network && events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE)) { + network_rebalance_necessary = true; + break; + } + } + network_rebalance_necessary = well_model_.comm().max(network_rebalance_necessary); + return network_rebalance_necessary; +} + +template +bool BlackoilWellModelNetworkGeneric:: +shouldBalance(const int reportStepIdx, const int iterationIdx) const +{ + // if network is not active, we do not need to balance the network + const auto& network = well_model_.schedule()[reportStepIdx].network(); + if (!network.active()) { + return false; + } + + const auto& balance = well_model_.schedule()[reportStepIdx].network_balance(); + if (balance.mode() == Network::Balance::CalcMode::TimeStepStart) { + return iterationIdx == 0; + } else if (balance.mode() == Network::Balance::CalcMode::NUPCOL) { + const int nupcol = well_model_.schedule()[reportStepIdx].nupcol(); + return iterationIdx < nupcol; + } else { + // We do not support any other rebalancing modes, + // i.e. TimeInterval based rebalancing is not available. + // This should be warned about elsewhere, so we choose to + // avoid spamming with a warning here. + return false; + } +} + +template +Scalar +BlackoilWellModelNetworkGeneric:: +updatePressures(const int reportStepIdx, + const Scalar damping_factor, + const Scalar upper_update_bound) +{ + OPM_TIMEFUNCTION(); + // Get the network and return if inactive (no wells in network at this time) + const auto& network = well_model_.schedule()[reportStepIdx].network(); + if (!network.active()) { + return 0.0; + } + + const auto previous_node_pressures = node_pressures_; + + node_pressures_ = this->computePressures(network, + *well_model_.getVFPProperties().getProd(), + reportStepIdx, + well_model_.comm()); + + // here, the network imbalance is the difference between the previous nodal pressure and the new nodal pressure + Scalar network_imbalance = 0.; + if (!this->active()) { + return network_imbalance; + } + + if (!previous_node_pressures.empty()) { + for (const auto& [name, new_pressure]: node_pressures_) { + if (previous_node_pressures.count(name) <= 0) { + if (std::abs(new_pressure) > network_imbalance) { + network_imbalance = std::abs(new_pressure); + } + continue; + } + const auto pressure = previous_node_pressures.at(name); + const Scalar change = (new_pressure - pressure); + if (std::abs(change) > network_imbalance) { + network_imbalance = std::abs(change); + } + // We dampen the nodal pressure change during one iteration since our nodal pressure calculation + // is somewhat explicit. There is a relative dampening factor applied to the update value, and also + // the maximum update is limited (to 5 bar by default, can be changed with --network-max-pressure-update-in-bars). + const Scalar damped_change = std::min(damping_factor * std::abs(change), upper_update_bound); + const Scalar sign = change > 0 ? 1. : -1.; + node_pressures_[name] = pressure + sign * damped_change; + } + } else { + for (const auto& [name, pressure]: node_pressures_) { + if (std::abs(pressure) > network_imbalance) { + network_imbalance = std::abs(pressure); + } + } + } + + for (auto& well : well_model_.genericWells()) { + + // Producers only, since we so far only support the + // "extended" network model (properties defined by + // BRANPROP and NODEPROP) which only applies to producers. + if (well->isProducer() && well->wellEcl().predictionMode()) { + const auto it = node_pressures_.find(well->wellEcl().groupName()); + if (it != node_pressures_.end()) { + // The well belongs to a group with has a network pressure constraint, + // set the dynamic THP constraint of the well accordingly. + const Scalar new_limit = it->second; + well->setDynamicThpLimit(new_limit); + SingleWellState& ws = well_model_.wellState()[well->indexOfWell()]; + const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP; + // TODO: not sure why the thp is NOT updated properly elsewhere + if (thp_is_limit) { + ws.thp = well->getTHPConstraint(well_model_.summaryState()); + } + } + } + } + return network_imbalance; +} + +template +void BlackoilWellModelNetworkGeneric:: +assignNodeValues(std::map& nodevalues, + const int reportStepIdx) const +{ + nodevalues.clear(); + if (reportStepIdx < 0) return; + + for (const auto& [node, pressure] : node_pressures_) { + nodevalues.emplace(node, data::NodeData{pressure}); + // Assign node values of well groups to GPR:WELLNAME + const auto& sched = well_model_.schedule(); + if (!sched.hasGroup(node, reportStepIdx)) { + continue; + } + const auto& group = sched.getGroup(node, reportStepIdx); + for (const std::string& wellname : group.wells()) { + nodevalues.emplace(wellname, data::NodeData{pressure}); + } + } + + const auto& network = well_model_.schedule()[reportStepIdx].network(); + if (!network.active()) { + return; + } + + auto converged_pressures = this->computePressures(network, + *well_model_.getVFPProperties().getProd(), + reportStepIdx, + well_model_.comm()); + for (const auto& [node, converged_pressure] : converged_pressures) { + auto it = nodevalues.find(node); + assert(it != nodevalues.end() ); + it->second.converged_pressure = converged_pressure; + // Assign node values of group to GPR:WELLNAME + const auto& sched = well_model_.schedule(); + if (!sched.hasGroup(node, reportStepIdx)) { + continue; + } + const auto& group = sched.getGroup(node, reportStepIdx); + for (const std::string& wellname : group.wells()) { + auto it2 = nodevalues.find(wellname); + assert(it2 != nodevalues.end()); + it2->second.converged_pressure = converged_pressure; + } + } +} + +template +void BlackoilWellModelNetworkGeneric:: +initialize(const int report_step) +{ + const auto& network = well_model_.schedule()[report_step].network(); + if (network.active() && !node_pressures_.empty()) { + for (auto& well : well_model_.genericWells()) { + initializeWell(*well); + } + } +} + +template +void BlackoilWellModelNetworkGeneric:: +initializeWell(WellInterfaceGeneric& well) +{ + // Producers only, since we so far only support the + // "extended" network model (properties defined by + // BRANPROP and NODEPROP) which only applies to producers. + if (well.isProducer() && !node_pressures_.empty()) { + const auto it = this->node_pressures_.find(well.wellEcl().groupName()); + if (it != this->node_pressures_.end()) { + // The well belongs to a group which has a network nodal pressure, + // set the dynamic THP constraint based on the network nodal pressure + well.setDynamicThpLimit(it->second); + } + } +} + +template +std::map +BlackoilWellModelNetworkGeneric:: +computePressures(const Network::ExtNetwork& network, + const VFPProdProperties& vfp_prod_props, + const int reportStepIdx, + const Parallel::Communication& comm) const +{ + // TODO: Only dealing with production networks for now. + OPM_TIMEFUNCTION(); + if (!network.active()) { + return {}; + } + + const auto& group_state = well_model_.groupStateHelper().groupState(); + const auto& well_state = well_model_.groupStateHelper().wellState(); + + std::map node_pressures; + auto roots = network.roots(); + for (const auto& root : roots) { + // Fixed pressure nodes of the network are the roots of trees. + // Leaf nodes must correspond to groups in the group structure. + // Let us first find all leaf nodes of the network. We also + // create a vector of all nodes, ordered so that a child is + // always after its parent. + std::stack children; + std::set leaf_nodes; + std::vector root_to_child_nodes; + // children.push(network.root().name()); + children.push(root.get().name()); + while (!children.empty()) { + const auto node = children.top(); + children.pop(); + root_to_child_nodes.push_back(node); + auto branches = network.downtree_branches(node); + if (branches.empty()) { + leaf_nodes.insert(node); + } + for (const auto& branch : branches) { + children.push(branch.downtree_node()); + } + } + assert(children.empty()); + + // Starting with the leaf nodes of the network, get the flow rates + // from the corresponding groups. + std::map> node_inflows; + const std::vector zero_rates(3, 0.0); + for (const auto& node : leaf_nodes) { + // Guard against empty leaf nodes (may not be present in GRUPTREE) + if (!well_model_.groupStateHelper().groupState().has_production_rates(node)) { + node_inflows[node] = zero_rates; + continue; + } + + node_inflows[node] = group_state.network_leaf_node_production_rates(node); + // Add the ALQ amounts to the gas rates if requested. + if (network.node(node).add_gas_lift_gas()) { + const auto& group = well_model_.schedule().getGroup(node, reportStepIdx); + Scalar alq = 0.0; + for (const std::string& wellname : group.wells()) { + const Well& well = well_model_.schedule().getWell(wellname, reportStepIdx); + if (well.isInjector() || !well_state.isOpen(wellname)) { + continue; + } + const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) + * well_state.getGlobalEfficiencyScalingFactor(wellname); + if (const auto& well_index = well_state.index(wellname); + well_index.has_value() && well_state.wellIsOwned(well_index.value(), wellname)) + { + alq += well_state.well(wellname).alq_state.get() * efficiency; + } + } + alq = comm.sum(alq); + node_inflows[node][IndexTraits::gasPhaseIdx] += alq; + } + } + + // Accumulate in the network, towards the roots. Note that a + // root (i.e. fixed pressure node) can still be contributing + // flow towards other nodes in the network, i.e. a node is + // the root of a subtree. + auto child_to_root_nodes = root_to_child_nodes; + std::reverse(child_to_root_nodes.begin(), child_to_root_nodes.end()); + for (const auto& node : child_to_root_nodes) { + const auto upbranch = network.uptree_branch(node); + if (upbranch) { + // Add downbranch rates to upbranch. + std::vector& up = node_inflows[(*upbranch).uptree_node()]; + const std::vector& down = node_inflows[node]; + // We now also support NEFAC + const Scalar efficiency = network.node(node).efficiency(); + if (up.empty()) { + up = std::vector(down.size(), 0.0); + } + assert(up.size() == down.size()); + for (std::size_t ii = 0; ii < up.size(); ++ii) { + up[ii] += efficiency * down[ii]; + } + } + } + + // Going the other way (from roots to leafs), calculate the pressure + // at each node using VFP tables and rates. + // std::map node_pressures; + for (const auto& node : root_to_child_nodes) { + auto press = network.node(node).terminal_pressure(); + if (press) { + node_pressures[node] = *press; + } + else { + const auto upbranch = network.uptree_branch(node); + assert(upbranch); + const Scalar up_press = node_pressures[(*upbranch).uptree_node()]; + const auto vfp_table = (*upbranch).vfp_table(); + if (vfp_table) { + OPM_TIMEBLOCK(NetworkVfpCalculations); + // The rates are here positive, but the VFP code expects the + // convention that production rates are negative, so we must + // take a copy and flip signs. + auto rates = node_inflows[node]; + std::transform( + rates.begin(), rates.end(), rates.begin(), [](const auto r) { return -r; }); + assert(rates.size() == 3); + // NB! ALQ in extended network is never implicitly the gas lift rate (GRAT), i.e., the + // gas lift rates only enters the network pressure calculations through the rates + // (e.g., in GOR calculations) unless a branch ALQ is set in BRANPROP. + // + // @TODO: Standard network + Scalar alq = (*upbranch).alq_value().value_or(0.0); + node_pressures[node] + = vfp_prod_props.bhp(*vfp_table, + rates[IndexTraits::waterPhaseIdx], + rates[IndexTraits::oilPhaseIdx], + rates[IndexTraits::gasPhaseIdx], + up_press, + alq, + 0.0, // explicit_wfr + 0.0, // explicit_gfr + false); // use_expvfp we dont support explicit lookup +#define EXTRA_DEBUG_NETWORK 0 +#if EXTRA_DEBUG_NETWORK + std::ostringstream oss; + oss << "parent: " << (*upbranch).uptree_node() << " child: " << node << " rates = [ " + << rates[0] * 86400 << ", " << rates[1] * 86400 << ", " << rates[2] * 86400 << " ]" + << " p(parent) = " << up_press / 1e5 << " p(child) = " << node_pressures[node] / 1e5 + << std::endl; + OpmLog::debug(oss.str()); +#endif + } + else { + // Table number specified as 9999 in the deck, no pressure loss. + if (network.node(node).as_choke()) { + // Node pressure is set to the group THP. + node_pressures[node] = well_model_.groupStateHelper().groupState().well_group_thp(node); + } + else { + node_pressures[node] = up_press; + } + } + } + } + } + return node_pressures; +} + +template +bool BlackoilWellModelNetworkGeneric:: +operator==(const BlackoilWellModelNetworkGeneric& rhs) const +{ + return + this->node_pressures_ == rhs.node_pressures_ + && this->last_valid_node_pressures_ == rhs.last_valid_node_pressures_; +} + +template class BlackoilWellModelNetworkGeneric; + +#if FLOW_INSTANTIATE_FLOAT +template class BlackoilWellModelNetworkGeneric; +#endif + +} diff --git a/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp new file mode 100644 index 00000000000..95989a6d259 --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelNetworkGeneric.hpp @@ -0,0 +1,125 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 - 2017 Statoil ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 IRIS AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BLACKOILWELLMODEL_NETWORK_GENERIC_HEADER_INCLUDED +#define OPM_BLACKOILWELLMODEL_NETWORK_GENERIC_HEADER_INCLUDED + +#include + +#include + +#include + +#include +#include +#include + +namespace Opm { + class Schedule; + template class BlackoilWellModelGeneric; + template class WellInterfaceGeneric; + template class VFPProdProperties; +} + +namespace Opm { + +/// Class for handling the blackoil well network model. +template +class BlackoilWellModelNetworkGeneric +{ +public: + BlackoilWellModelNetworkGeneric(BlackoilWellModelGeneric& well_model); + + virtual ~BlackoilWellModelNetworkGeneric() = default; + + /// return true if network is active (at least one network well in prediction mode) + bool active() const + { return active_; } + + const std::map& + nodePressures() const { return node_pressures_; } + + // do not use, only needed for serialization testing + void setNodePressures(const std::map& values) + { node_pressures_ = values; } + + void setFromRestart(const std::optional>& restart_pressures); + + //! \brief Initialize wells according to network configuration. + void initialize(const int report_step); + + //! \brief Initialize a single well according to network configuration. + void initializeWell(WellInterfaceGeneric& well); + + /// Checks if network is active (at least one network well on prediction). + void updateActiveState(const int report_step); + + /// Checks if there are reasons to perform a pre-step network re-balance. + /// (Currently, the only reasons are network well status changes.) + /// (TODO: Consider if adding network change events would be helpful.) + bool needPreStepRebalance(const int report_step) const; + + bool shouldBalance(const int reportStepIndex, + const int iterationIdx) const; + + Scalar updatePressures(const int reportStepIdx, + const Scalar damping_factor, + const Scalar update_upper_bound); + + void assignNodeValues(std::map& nodevalues, + const int reportStepIdx) const; + + void commitState() + { this->last_valid_node_pressures_ = this->node_pressures_; } + + void resetState() + { this->node_pressures_ = this->last_valid_node_pressures_; } + + template + void serializeOp(Serializer& serializer) + { + serializer(node_pressures_); + serializer(last_valid_node_pressures_); + } + + bool operator==(const BlackoilWellModelNetworkGeneric& rhs) const; + +protected: + std::map + computePressures(const Network::ExtNetwork& network, + const VFPProdProperties& vfp_prod_props, + const int reportStepIdx, + const Parallel::Communication& comm) const; + + + bool active_{false}; + BlackoilWellModelGeneric& well_model_; + + // Network pressures for output and initialization + std::map node_pressures_; + // Valid network pressures for output and initialization for safe restart after failed iterations + std::map last_valid_node_pressures_; +}; + +} // namespace Opm + +#endif diff --git a/opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp b/opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp new file mode 100644 index 00000000000..d9b6e446810 --- /dev/null +++ b/opm/simulators/wells/BlackoilWellModelNetwork_impl.hpp @@ -0,0 +1,384 @@ +/* + Copyright 2016 - 2019 SINTEF Digital, Mathematics & Cybernetics. + Copyright 2016 - 2018 Equinor ASA. + Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2016 - 2018 Norce AS + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_BLACKOILWELLMODEL_NETWORK_IMPL_HEADER_INCLUDED +#define OPM_BLACKOILWELLMODEL_NETWORK_IMPL_HEADER_INCLUDED + +// Improve IDE experience +#ifndef OPM_BLACKOILWELLMODEL_NETWORK_HEADER_INCLUDED +#include +#include +#endif + +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include + +namespace Opm { + +template +BlackoilWellModelNetwork:: +BlackoilWellModelNetwork(BlackoilWellModel& well_model) + : BaseType(well_model) + , well_model_(well_model) +{} + +template +void +BlackoilWellModelNetwork:: +doPreStepRebalance(DeferredLogger& deferred_logger) +{ + OPM_TIMEFUNCTION(); + const double dt = well_model_.simulator().timeStepSize(); + // TODO: should we also have the group and network backed-up + // here in case the solution did not get converged? + auto& well_state = well_model_.wellState(); + + const bool changed_well_group = + well_model_.updateWellControlsAndNetwork(true, dt, deferred_logger); + well_model_.assembleWellEqWithoutIteration(dt, deferred_logger); + const bool converged = + well_model_.getWellConvergence(well_model_.B_avg(), true).converged() && + !changed_well_group; + + OPM_BEGIN_PARALLEL_TRY_CATCH(); + for (auto& well : this->well_model_) { + well->solveEqAndUpdateWellState(well_model_.simulator(), + well_state, + deferred_logger); + } + OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModelNetwork::doPreStepRebalance() failed: ", + well_model_.simulator().vanguard().grid().comm()); + + if (!converged) { + const std::string msg = + fmt::format("Initial (pre-step) network balance did not converge."); + deferred_logger.warning(msg); + } +} + +template +std::tuple::Scalar> +BlackoilWellModelNetwork:: +update(const bool mandatory_network_balance, + DeferredLogger& deferred_logger, + const bool relax_network_tolerance) +{ + OPM_TIMEFUNCTION(); + const int episodeIdx = well_model_.simulator().episodeIndex(); + const auto& network = well_model_.schedule()[episodeIdx].network(); + if (!well_model_.wellsActive() && !network.active()) { + return {false, 0.0}; + } + + const int iterationIdx = well_model_.simulator().model().newtonMethod().numIterations(); + const auto& comm = well_model_.simulator().vanguard().grid().comm(); + + // network related + Scalar network_imbalance = 0.0; + bool more_network_update = false; + if (this->shouldBalance(episodeIdx, iterationIdx) || mandatory_network_balance) { + OPM_TIMEBLOCK(BalanceNetwork); + const double dt = well_model_.simulator().timeStepSize(); + // Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES) + const bool well_group_thp_updated = computeWellGroupThp(dt, deferred_logger); + const int max_number_of_sub_iterations = + well_model_.param().network_max_sub_iterations_; + const Scalar network_pressure_update_damping_factor = + well_model_.param().network_pressure_update_damping_factor_; + const Scalar network_max_pressure_update = + well_model_.param().network_max_pressure_update_in_bars_ * unit::barsa; + bool more_network_sub_update = false; + for (int i = 0; i < max_number_of_sub_iterations; i++) { + const auto local_network_imbalance = + this->updatePressures(episodeIdx, + network_pressure_update_damping_factor, + network_max_pressure_update); + network_imbalance = comm.max(local_network_imbalance); + const auto& balance = well_model_.schedule()[episodeIdx].network_balance(); + constexpr Scalar relaxation_factor = 10.0; + const Scalar tolerance = + relax_network_tolerance ? relaxation_factor * balance.pressure_tolerance() + : balance.pressure_tolerance(); + more_network_sub_update = this->active() && network_imbalance > tolerance; + if (!more_network_sub_update) { + break; + } + + for (const auto& well : well_model_) { + if (well->isInjector() || !well->wellEcl().predictionMode()) { + continue; + } + + const auto it = this->node_pressures_.find(well->wellEcl().groupName()); + if (it != this->node_pressures_.end()) { + well->prepareWellBeforeAssembling(well_model_.simulator(), + dt, + well_model_.groupStateHelper(), + well_model_.wellState(), + deferred_logger); + } + } + well_model_.updateAndCommunicateGroupData(episodeIdx, + iterationIdx, + well_model_.param().nupcol_group_rate_tolerance_, + /*update_wellgrouptarget*/ true, + deferred_logger); + } + more_network_update = more_network_sub_update || well_group_thp_updated; + } + return { more_network_update, network_imbalance }; +} + +template +bool +BlackoilWellModelNetwork:: +computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger) +{ + OPM_TIMEFUNCTION(); + const int reportStepIdx = well_model_.simulator().episodeIndex(); + const auto& network = well_model_.schedule()[reportStepIdx].network(); + const auto& balance = well_model_.schedule()[reportStepIdx].network_balance(); + const Scalar thp_tolerance = balance.thp_tolerance(); + + if (!network.active()) { + return false; + } + + auto& well_state = well_model_.wellState(); + auto& group_state = well_model_.groupState(); + + bool well_group_thp_updated = false; + for (const std::string& nodeName : network.node_names()) { + const bool has_choke = network.node(nodeName).as_choke(); + if (has_choke) { + const auto& summary_state = well_model_.simulator().vanguard().summaryState(); + const Group& group = well_model_.schedule().getGroup(nodeName, reportStepIdx); + + //TODO: Auto choke combined with RESV control is not supported + std::vector resv_coeff(Indices::numPhases, 1.0); + Scalar gratTargetFromSales = 0.0; + if (group_state.has_grat_sales_target(group.name())) + gratTargetFromSales = group_state.grat_sales_target(group.name()); + + const auto ctrl = group.productionControls(summary_state); + auto cmode_tmp = ctrl.cmode; + Scalar target_tmp{0.0}; + bool fld_none = false; + if (cmode_tmp == Group::ProductionCMode::FLD || cmode_tmp == Group::ProductionCMode::NONE) { + fld_none = true; + // Target is set for an ancestor group. Target for autochoke group to be + // derived via group guide rates + const Scalar efficiencyFactor = 1.0; + const Group& parentGroup = well_model_.schedule().getGroup(group.parent(), reportStepIdx); + auto target = WellGroupControls:: + getAutoChokeGroupProductionTargetRate(group.name(), + parentGroup, + well_model_.groupStateHelper(), + well_model_.schedule(), + summary_state, + resv_coeff, + efficiencyFactor, + reportStepIdx, + &well_model_.guideRate(), + local_deferredLogger); + target_tmp = target.first; + cmode_tmp = target.second; + } + const auto cmode = cmode_tmp; + using TargetCalculatorType = GroupStateHelpers::TargetCalculator; + TargetCalculatorType tcalc(cmode, FluidSystem::phaseUsage(), resv_coeff, + gratTargetFromSales, nodeName, group_state, + group.has_gpmaint_control(cmode)); + if (!fld_none) + { + // Target is set for the autochoke group itself + target_tmp = tcalc.groupTarget(ctrl, local_deferredLogger); + } + + const Scalar orig_target = target_tmp; + + auto mismatch = [&] (auto group_thp) { + Scalar group_rate(0.0); + Scalar rate(0.0); + for (auto& well : well_model_) { + std::string well_name = well->name(); + auto& ws = well_state.well(well_name); + if (group.hasWell(well_name)) { + well->setDynamicThpLimit(group_thp); + const Well& well_ecl = well_model_.eclWells()[well->indexOfWell()]; + const auto inj_controls = Well::InjectionControls(0); + const auto prod_controls = well_ecl.productionControls(summary_state); + well->iterateWellEqWithSwitching(well_model_.simulator(), + dt, + inj_controls, + prod_controls, + well_model_.groupStateHelper(), + well_state, + local_deferredLogger, + /*fixed_control=*/false, + /*fixed_status=*/false, + /*sovling_with_zero_rate=*/false); + rate = -tcalc.calcModeRateFromRates(ws.surface_rates); + group_rate += rate; + } + } + return (group_rate - orig_target)/orig_target; + }; + + 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); + } + + using WellBhpThpCalculatorType = WellBhpThpCalculator; + //Find an initial bracket + 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(); + WellBhpThpCalculatorType::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp); + // Narrow down the bracket + Scalar low1, high1; + std::array range = {Scalar{0.9}*min_thp, Scalar{1.1}*max_thp}; + std::optional appr_sol; + WellBhpThpCalculatorType::bruteForceBracketCommonTHP(mismatch, + range, + low1, + high1, + appr_sol, + 0.0, + local_deferredLogger); + min_thp = low1; + max_thp = high1; + range_initial = {min_thp, max_thp}; + } + + 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{Scalar{0.9} * autochoke_thp.value(), + Scalar{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 group THP"); + const bool finding_bracket = WellBhpThpCalculatorType:: + bruteForceBracketCommonTHP(mismatch, + range, + low, + high, + approximate_solution, + tolerance1, + local_deferredLogger); + + if (approximate_solution.has_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; + 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("Group THP value = " + std::to_string(autochoke_thp.value())); + } else { + autochoke_thp.reset(); + local_deferredLogger.debug("Group THP solve failed due to bracketing failure"); + } + } + 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 : well_model_) { + std::string well_name = well->name(); + + if (well->isInjector() || !well->wellEcl().predictionMode()) + continue; + + if (group.hasWell(well_name)) { + well->setDynamicThpLimit(well_group_thp); + } + const auto& ws = well_model_.wellState().well(well->indexOfWell()); + const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP; + if (thp_is_limit) { + well->prepareWellBeforeAssembling(well_model_.simulator(), + dt, + well_model_.groupStateHelper(), + well_model_.wellState(), + local_deferredLogger); + } + } + + // Use the group THP in computeNetworkPressures(). + const auto& current_well_group_thp = group_state.is_autochoke_group(nodeName) + ? group_state.well_group_thp(nodeName) + : 1e30; + if (std::abs(current_well_group_thp - well_group_thp) > balance.pressure_tolerance()) { + well_group_thp_updated = true; + group_state.update_well_group_thp(nodeName, well_group_thp); + } + } + } + return well_group_thp_updated; +} + +} // namespace Opm + +#endif // OPM_BLACKOILWELLMODEL_NETWORK_IMPL_HEADER_INCLUDED diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 30bd08e4da9..b1e6ca65db4 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -30,7 +30,6 @@ #endif #include -#include #include #include @@ -47,10 +46,7 @@ #include #include #include -#include -#include #include -#include #ifdef RESERVOIR_COUPLING_ENABLED #include @@ -78,11 +74,12 @@ namespace Opm { BlackoilWellModel(Simulator& simulator) : WellConnectionModule(*this, simulator.gridView().comm()) , BlackoilWellModelGeneric(simulator.vanguard().schedule(), - gaslift_, - simulator.vanguard().summaryState(), - simulator.vanguard().eclState(), - FluidSystem::phaseUsage(), - simulator.gridView().comm()) + gaslift_, + network_, + simulator.vanguard().summaryState(), + simulator.vanguard().eclState(), + FluidSystem::phaseUsage(), + simulator.gridView().comm()) , simulator_(simulator) , guide_rate_handler_{ *this, @@ -91,6 +88,7 @@ namespace Opm { simulator.vanguard().grid().comm() } , gaslift_(this->terminal_output_) + , network_(*this) { local_num_cells_ = simulator_.gridView().size(0); @@ -626,16 +624,8 @@ namespace Opm { } const auto& network = this->schedule()[timeStepIdx].network(); - if (network.active() && !this->node_pressures_.empty()) { - if (well->isProducer()) { - const auto it = this->node_pressures_.find(well->wellEcl().groupName()); - if (it != this->node_pressures_.end()) { - // The well belongs to a group which has a network nodal pressure, - // set the dynamic THP constraint based on the network nodal pressure - const Scalar nodal_pressure = it->second; - well->setDynamicThpLimit(nodal_pressure); - } - } + if (network.active()) { + this->network_.initializeWell(*well); } try { using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells; @@ -1048,27 +1038,12 @@ namespace Opm { } this->well_container_generic_.clear(); - for (auto& w : well_container_) + for (auto& w : well_container_) { this->well_container_generic_.push_back(w.get()); - - const auto& network = this->schedule()[report_step].network(); - if (network.active() && !this->node_pressures_.empty()) { - for (auto& well: this->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. - if (well->isProducer()) { - const auto it = this->node_pressures_.find(well->wellEcl().groupName()); - if (it != this->node_pressures_.end()) { - // The well belongs to a group which has a network nodal pressure, - // set the dynamic THP constraint based on the network nodal pressure - const Scalar nodal_pressure = it->second; - well->setDynamicThpLimit(nodal_pressure); - } - } - } } + this->network_.initialize(report_step); + this->wbp_.registerOpenWellsForWBPCalculation(); } @@ -1152,36 +1127,6 @@ namespace Opm { - template - void - BlackoilWellModel:: - doPreStepNetworkRebalance(DeferredLogger& deferred_logger) - { - OPM_TIMEFUNCTION(); - const double dt = this->simulator_.timeStepSize(); - // TODO: should we also have the group and network backed-up here in case the solution did not get converged? - auto& well_state = this->wellState(); - - const bool changed_well_group = updateWellControlsAndNetwork(true, dt, deferred_logger); - assembleWellEqWithoutIteration(dt, deferred_logger); - const bool converged = this->getWellConvergence(this->B_avg_, true).converged() && !changed_well_group; - - OPM_BEGIN_PARALLEL_TRY_CATCH(); - for (auto& well : this->well_container_) { - well->solveEqAndUpdateWellState(simulator_, well_state, deferred_logger); - } - OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::doPreStepNetworkRebalance() failed: ", - this->simulator_.vanguard().grid().comm()); - - if (!converged) { - const std::string msg = fmt::format("Initial (pre-step) network balance did not converge."); - deferred_logger.warning(msg); - } - } - - - - template void BlackoilWellModel:: @@ -1272,7 +1217,7 @@ namespace Opm { // only output to terminal if we at the last newton iterations where we try to balance the network. const int episodeIdx = simulator_.episodeIndex(); const int iterationIdx = simulator_.model().newtonMethod().numIterations(); - if (this->shouldBalanceNetwork(episodeIdx, iterationIdx + 1)) { + if (this->network_.shouldBalance(episodeIdx, iterationIdx + 1)) { if (this->terminal_output_) { const std::string msg = fmt::format("Maximum of {:d} network iterations has been used and we stop the update, \n" "and try again after the next Newton iteration (imbalance = {:.2e} bar, ctrl_change = {})", @@ -1328,10 +1273,9 @@ namespace Opm { // and in prepareWellsBeforeAssembling during well solves. bool well_group_control_changed = updateWellControls(local_deferredLogger); const auto [more_inner_network_update, network_imbalance] = - updateNetworks(mandatory_network_balance, - local_deferredLogger, - relax_network_tolerance); - + this->network_.update(mandatory_network_balance, + local_deferredLogger, + relax_network_tolerance); bool alq_updated = false; OPM_BEGIN_PARALLEL_TRY_CATCH(); @@ -1339,10 +1283,11 @@ namespace Opm { if (optimize_gas_lift) { // we need to update the potentials if the thp limit as been modified by // the network balancing - const bool updatePotentials = (this->shouldBalanceNetwork(reportStepIdx, iterationIdx) || mandatory_network_balance); + const bool updatePotentials = (this->network_.shouldBalance(reportStepIdx, iterationIdx) || + mandatory_network_balance); alq_updated = gaslift_.maybeDoGasLiftOptimize(simulator_, well_container_, - this->node_pressures_, + this->network_.nodePressures(), updatePotentials, this->wellState(), this->groupState(), @@ -1367,204 +1312,11 @@ namespace Opm { } // we need to re-iterate the network when the well group controls changed or gaslift/alq is changed or // the inner iterations are did not converge - const bool more_network_update = this->shouldBalanceNetwork(reportStepIdx, iterationIdx) && + const bool more_network_update = this->network_.shouldBalance(reportStepIdx, iterationIdx) && (more_inner_network_update || well_group_control_changed || alq_updated); return {well_group_control_changed, more_network_update, network_imbalance}; } - // 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 - bool - BlackoilWellModel:: - computeWellGroupThp(const double dt, DeferredLogger& local_deferredLogger) - { - OPM_TIMEFUNCTION(); - 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 false; - } - - auto& well_state = this->wellState(); - auto& group_state = this->groupState(); - - bool well_group_thp_updated = false; - for (const std::string& nodeName : network.node_names()) { - const bool has_choke = network.node(nodeName).as_choke(); - if (has_choke) { - const auto& summary_state = this->simulator_.vanguard().summaryState(); - const Group& group = this->schedule().getGroup(nodeName, reportStepIdx); - - //TODO: Auto choke combined with RESV control is not supported - std::vector resv_coeff(Indices::numPhases, 1.0); - Scalar gratTargetFromSales = 0.0; - if (group_state.has_grat_sales_target(group.name())) - gratTargetFromSales = group_state.grat_sales_target(group.name()); - - const auto ctrl = group.productionControls(summary_state); - auto cmode_tmp = ctrl.cmode; - Scalar target_tmp{0.0}; - bool fld_none = false; - if (cmode_tmp == Group::ProductionCMode::FLD || cmode_tmp == Group::ProductionCMode::NONE) { - fld_none = true; - // Target is set for an ancestor group. Target for autochoke group to be - // derived via group guide rates - const Scalar efficiencyFactor = 1.0; - const Group& parentGroup = this->schedule().getGroup(group.parent(), reportStepIdx); - auto target = WellGroupControls::getAutoChokeGroupProductionTargetRate( - group.name(), - parentGroup, - this->groupStateHelper(), - this->schedule(), - summary_state, - resv_coeff, - efficiencyFactor, - reportStepIdx, - &this->guideRate_, - local_deferredLogger); - target_tmp = target.first; - cmode_tmp = target.second; - } - const auto cmode = cmode_tmp; - using TargetCalculatorType = GroupStateHelpers -::TargetCalculator; - TargetCalculatorType tcalc(cmode, FluidSystem::phaseUsage(), resv_coeff, - gratTargetFromSales, nodeName, group_state, - group.has_gpmaint_control(cmode)); - if (!fld_none) - { - // Target is set for the autochoke group itself - target_tmp = tcalc.groupTarget(ctrl, local_deferredLogger); - } - - const Scalar orig_target = target_tmp; - - 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(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, this->groupStateHelper(), this->wellState(), local_deferredLogger, - /*fixed_control=*/false, - /*fixed_status=*/false, - /*solving_with_zero_rate=*/false - ); - rate = -tcalc.calcModeRateFromRates(ws.surface_rates); - group_rate += rate; - } - } - return (group_rate - orig_target)/orig_target; - }; - - 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); - } - - using WellBhpThpCalculatorType = WellBhpThpCalculator; - //Find an initial bracket - 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(); - WellBhpThpCalculatorType::bruteForceBracketCommonTHP(mismatch, min_thp, max_thp); - // Narrow down the bracket - Scalar low1, high1; - std::array range = {Scalar{0.9}*min_thp, Scalar{1.1}*max_thp}; - std::optional appr_sol; - WellBhpThpCalculatorType::bruteForceBracketCommonTHP(mismatch, range, low1, high1, appr_sol, 0.0, local_deferredLogger); - min_thp = low1; - max_thp = high1; - range_initial = {min_thp, max_thp}; - } - - 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{Scalar{0.9} * autochoke_thp.value(), - Scalar{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 group THP"); - const bool finding_bracket = WellBhpThpCalculatorType::bruteForceBracketCommonTHP(mismatch, range, low, high, approximate_solution, tolerance1, local_deferredLogger); - - if (approximate_solution.has_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; - 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("Group THP value = " + std::to_string(autochoke_thp.value())); - } else { - autochoke_thp.reset(); - local_deferredLogger.debug("Group THP solve failed due to bracketing failure"); - } - } - 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 (well->isInjector() || !well->wellEcl().predictionMode()) - continue; - - if (group.hasWell(well_name)) { - well->setDynamicThpLimit(well_group_thp); - } - const auto& ws = this->wellState().well(well->indexOfWell()); - const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP; - if (thp_is_limit) { - well->prepareWellBeforeAssembling( - this->simulator_, dt, this->groupStateHelper(), this->wellState(), local_deferredLogger - ); - } - } - - // Use the group THP in computeNetworkPressures(). - const auto& current_well_group_thp = group_state.is_autochoke_group(nodeName) ? group_state.well_group_thp(nodeName) : 1e30; - if (std::abs(current_well_group_thp - well_group_thp) > balance.pressure_tolerance()) { - well_group_thp_updated = true; - group_state.update_well_group_thp(nodeName, well_group_thp); - } - } - } - return well_group_thp_updated; - } - template void BlackoilWellModel:: @@ -1947,63 +1699,6 @@ ::TargetCalculator; } - template - std::tuple::Scalar> - BlackoilWellModel:: - updateNetworks(const bool mandatory_network_balance, - DeferredLogger& deferred_logger, - const bool relax_network_tolerance) - { - OPM_TIMEFUNCTION(); - const int episodeIdx = simulator_.episodeIndex(); - const auto& network = this->schedule()[episodeIdx].network(); - if (!this->wellsActive() && !network.active()) { - return {false, 0.0}; - } - - const int iterationIdx = simulator_.model().newtonMethod().numIterations(); - const auto& comm = simulator_.vanguard().grid().comm(); - - // network related - Scalar network_imbalance = 0.0; - bool more_network_update = false; - if (this->shouldBalanceNetwork(episodeIdx, iterationIdx) || mandatory_network_balance) { - OPM_TIMEBLOCK(BalanceNetwork); - const double dt = this->simulator_.timeStepSize(); - // Calculate common THP for subsea manifold well group (item 3 of NODEPROP set to YES) - const bool well_group_thp_updated = computeWellGroupThp(dt, deferred_logger); - const int max_number_of_sub_iterations = param_.network_max_sub_iterations_; - const Scalar network_pressure_update_damping_factor = param_.network_pressure_update_damping_factor_; - const Scalar network_max_pressure_update = param_.network_max_pressure_update_in_bars_ * unit::barsa; - bool more_network_sub_update = false; - for (int i = 0; i < max_number_of_sub_iterations; i++) { - const auto local_network_imbalance = this->updateNetworkPressures(episodeIdx, network_pressure_update_damping_factor, network_max_pressure_update); - network_imbalance = comm.max(local_network_imbalance); - const auto& balance = this->schedule()[episodeIdx].network_balance(); - constexpr Scalar relaxation_factor = 10.0; - const Scalar tolerance = relax_network_tolerance ? relaxation_factor * balance.pressure_tolerance() : balance.pressure_tolerance(); - more_network_sub_update = this->networkActive() && network_imbalance > tolerance; - if (!more_network_sub_update) - break; - - for (const auto& well : well_container_) { - if (well->isInjector() || !well->wellEcl().predictionMode()) - continue; - - const auto it = this->node_pressures_.find(well->wellEcl().groupName()); - if (it != this->node_pressures_.end()) { - well->prepareWellBeforeAssembling(this->simulator_, dt, this->groupStateHelper(), this->wellState(), deferred_logger); - } - } - this->updateAndCommunicateGroupData(episodeIdx, iterationIdx, param_.nupcol_group_rate_tolerance_, - /*update_wellgrouptarget*/ true, deferred_logger); - } - more_network_update = more_network_sub_update || well_group_thp_updated; - } - return { more_network_update, network_imbalance }; - } - - template void BlackoilWellModel:: @@ -2268,11 +1963,12 @@ ::TargetCalculator; { // Check if there is a network with active prediction wells at this time step. const auto episodeIdx = simulator_.episodeIndex(); - this->updateNetworkActiveState(episodeIdx); + this->network_.updateActiveState(episodeIdx); // Rebalance the network initially if any wells in the network have status changes // (Need to check this before clearing events) - const bool do_prestep_network_rebalance = param_.pre_solve_network_ && this->needPreStepNetworkRebalance(episodeIdx); + const bool do_prestep_network_rebalance = + param_.pre_solve_network_ && this->network_.needPreStepRebalance(episodeIdx); for (const auto& well : well_container_) { auto& events = this->wellState().well(well->indexOfWell()).events; @@ -2307,7 +2003,9 @@ ::TargetCalculator; updatePrimaryVariables(deferred_logger); // Actually do the pre-step network rebalance, using the updated well states and initial solutions - if (do_prestep_network_rebalance) doPreStepNetworkRebalance(deferred_logger); + if (do_prestep_network_rebalance) { + network_.doPreStepRebalance(deferred_logger); + } } template diff --git a/opm/simulators/wells/GroupStateHelper.cpp b/opm/simulators/wells/GroupStateHelper.cpp index 93f0b77c4f6..cddb2cc15f4 100644 --- a/opm/simulators/wells/GroupStateHelper.cpp +++ b/opm/simulators/wells/GroupStateHelper.cpp @@ -419,163 +419,6 @@ GroupStateHelper::checkGroupConstraintsProd(const std::stri return std::make_pair(current_rate_available > target_rate_available, scale); } -template -std::map -GroupStateHelper::computeNetworkPressures(const Network::ExtNetwork& network, - const VFPProdProperties& vfp_prod_props, - const Parallel::Communication& comm) const -{ - // TODO: Only dealing with production networks for now. - OPM_TIMEFUNCTION(); - if (!network.active()) { - return {}; - } - - std::map node_pressures; - auto roots = network.roots(); - for (const auto& root : roots) { - // Fixed pressure nodes of the network are the roots of trees. - // Leaf nodes must correspond to groups in the group structure. - // Let us first find all leaf nodes of the network. We also - // create a vector of all nodes, ordered so that a child is - // always after its parent. - std::stack children; - std::set leaf_nodes; - std::vector root_to_child_nodes; - // children.push(network.root().name()); - children.push(root.get().name()); - while (!children.empty()) { - const auto node = children.top(); - children.pop(); - root_to_child_nodes.push_back(node); - auto branches = network.downtree_branches(node); - if (branches.empty()) { - leaf_nodes.insert(node); - } - for (const auto& branch : branches) { - children.push(branch.downtree_node()); - } - } - assert(children.empty()); - - // Starting with the leaf nodes of the network, get the flow rates - // from the corresponding groups. - std::map> node_inflows; - const std::vector zero_rates(3, 0.0); - for (const auto& node : leaf_nodes) { - // Guard against empty leaf nodes (may not be present in GRUPTREE) - if (!this->groupState().has_production_rates(node)) { - node_inflows[node] = zero_rates; - continue; - } - - node_inflows[node] = this->groupState().network_leaf_node_production_rates(node); - // Add the ALQ amounts to the gas rates if requested. - if (network.node(node).add_gas_lift_gas()) { - const auto& group = this->schedule_.getGroup(node, this->report_step_); - Scalar alq = 0.0; - for (const std::string& wellname : group.wells()) { - const Well& well = this->schedule_.getWell(wellname, this->report_step_); - if (well.isInjector() || !this->wellState().isOpen(wellname)) - continue; - const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) - * this->wellState().getGlobalEfficiencyScalingFactor(wellname); - const auto& well_index = this->wellState().index(wellname); - if (well_index.has_value() - && this->wellState().wellIsOwned(well_index.value(), wellname)) { - alq += this->wellState().well(wellname).alq_state.get() * efficiency; - } - } - alq = comm.sum(alq); - node_inflows[node][IndexTraits::gasPhaseIdx] += alq; - } - } - - // Accumulate in the network, towards the roots. Note that a - // root (i.e. fixed pressure node) can still be contributing - // flow towards other nodes in the network, i.e. a node is - // the root of a subtree. - auto child_to_root_nodes = root_to_child_nodes; - std::reverse(child_to_root_nodes.begin(), child_to_root_nodes.end()); - for (const auto& node : child_to_root_nodes) { - const auto upbranch = network.uptree_branch(node); - if (upbranch) { - // Add downbranch rates to upbranch. - std::vector& up = node_inflows[(*upbranch).uptree_node()]; - const std::vector& down = node_inflows[node]; - // We now also support NEFAC - const Scalar efficiency = network.node(node).efficiency(); - if (up.empty()) { - up = std::vector(down.size(), 0.0); - } - assert(up.size() == down.size()); - for (std::size_t ii = 0; ii < up.size(); ++ii) { - up[ii] += efficiency * down[ii]; - } - } - } - - // Going the other way (from roots to leafs), calculate the pressure - // at each node using VFP tables and rates. - // std::map node_pressures; - for (const auto& node : root_to_child_nodes) { - auto press = network.node(node).terminal_pressure(); - if (press) { - node_pressures[node] = *press; - } else { - const auto upbranch = network.uptree_branch(node); - assert(upbranch); - const Scalar up_press = node_pressures[(*upbranch).uptree_node()]; - const auto vfp_table = (*upbranch).vfp_table(); - if (vfp_table) { - OPM_TIMEBLOCK(NetworkVfpCalculations); - // The rates are here positive, but the VFP code expects the - // convention that production rates are negative, so we must - // take a copy and flip signs. - auto rates = node_inflows[node]; - std::transform( - rates.begin(), rates.end(), rates.begin(), [](const auto r) { return -r; }); - assert(rates.size() == 3); - // NB! ALQ in extended network is never implicitly the gas lift rate (GRAT), i.e., the - // gas lift rates only enters the network pressure calculations through the rates - // (e.g., in GOR calculations) unless a branch ALQ is set in BRANPROP. - // - // @TODO: Standard network - Scalar alq = (*upbranch).alq_value().value_or(0.0); - node_pressures[node] - = vfp_prod_props.bhp(*vfp_table, - rates[IndexTraits::waterPhaseIdx], - rates[IndexTraits::oilPhaseIdx], - rates[IndexTraits::gasPhaseIdx], - up_press, - alq, - 0.0, // explicit_wfr - 0.0, // explicit_gfr - false); // use_expvfp we dont support explicit lookup -#define EXTRA_DEBUG_NETWORK 0 -#if EXTRA_DEBUG_NETWORK - std::ostringstream oss; - oss << "parent: " << (*upbranch).uptree_node() << " child: " << node << " rates = [ " - << rates[0] * 86400 << ", " << rates[1] * 86400 << ", " << rates[2] * 86400 << " ]" - << " p(parent) = " << up_press / 1e5 << " p(child) = " << node_pressures[node] / 1e5 - << std::endl; - OpmLog::debug(oss.str()); -#endif - } else { - // Table number specified as 9999 in the deck, no pressure loss. - if (network.node(node).as_choke()) { - // Node pressure is set to the group THP. - node_pressures[node] = this->groupState().well_group_thp(node); - } else { - node_pressures[node] = up_press; - } - } - } - } - } - return node_pressures; -} - template Scalar GroupStateHelper::getGuideRate(const std::string& name, diff --git a/opm/simulators/wells/GroupStateHelper.hpp b/opm/simulators/wells/GroupStateHelper.hpp index 6348c50988a..d9b16645882 100644 --- a/opm/simulators/wells/GroupStateHelper.hpp +++ b/opm/simulators/wells/GroupStateHelper.hpp @@ -144,10 +144,6 @@ class GroupStateHelper const bool check_guide_rate, DeferredLogger& deferred_logger) const; - std::map computeNetworkPressures(const Network::ExtNetwork& network, - const VFPProdProperties& vfp_prod_props, - const Parallel::Communication& comm) const; - Scalar getGuideRate(const std::string& name, const GuideRateModel::Target target) const; GuideRate::RateVector getProductionGroupRateVector(const std::string& group_name) const; diff --git a/tests/test_RestartSerialization.cpp b/tests/test_RestartSerialization.cpp index 2ffff9198ff..20616fafa73 100644 --- a/tests/test_RestartSerialization.cpp +++ b/tests/test_RestartSerialization.cpp @@ -310,8 +310,9 @@ class BlackoilWellModelGenericTest : public BlackoilWellModelGeneric(schedule, gaslift, summaryState, + : BlackoilWellModelGeneric(schedule, gaslift, network_, summaryState, eclState, phase_usage, comm) + , network_(*this) { if (deserialize) { active_wgstate_.well_state = WellState(dummy); @@ -328,7 +329,7 @@ class BlackoilWellModelGenericTest : public BlackoilWellModelGeneric::serializationTestObject(dummy); last_valid_wgstate_ = WGState::serializationTestObject(dummy); nupcol_wgstate_ = WGState::serializationTestObject(dummy); @@ -371,6 +372,7 @@ class BlackoilWellModelGenericTest : public BlackoilWellModelGeneric network_; ParallelWellInfo dummy; };