diff --git a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp index bd54deaa0..32caff886 100644 --- a/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp +++ b/opm/simulators/utils/PartiallySupportedFlowKeywords.cpp @@ -381,6 +381,12 @@ partiallySupported() {13,{false, allow_values {10}, "MESSAGES(PRTGRPMS): option is not supported – will continue"}}, // GROUP_PRINT_LIMIT }, }, + { + "NETBALAN", + { + {5,{false, allow_values {10}, "NETBALAN(THPMXITE): option is not supported – will continue"}}, // MAX_ITER_THP + }, + }, { "NETWORK", { @@ -579,6 +585,16 @@ partiallySupported() {3,{false, allow_values {}, "MULTFLT(FLT-DIF): the diffusivity multiplier option is not supported – will continue"}}, // NOT_DEFINED }, }, + { + "NETBALAN", + { + {1,{true, [](const double value) { return value <= 0.0; }, "NETBALAN(NSTEP): only negative values or 0 supported – will stop"}}, // TIME_INTERVAL + {4,{false, allow_values {0.01}, "NETBALAN(GRPCNV): not supported – will continue"}}, // THP_CONVERGENCE_LIMIT + {6,{false, allow_values {1e20}, "NETBALAN(NTRGERR): not supported – will continue"}}, // TARGET_BALANCE_ERROR + {7,{false, allow_values {1e20}, "NETBALAN(NMAXERR): not supported – will continue"}}, // MAX_BALANCE_ERROR + {8,{false, allow_values {}, "NETBALAN(NTSMIN): not supported – will continue"}}, // MIN_TIME_STEP + }, + }, { "NNC", { diff --git a/opm/simulators/utils/UnsupportedFlowKeywords.cpp b/opm/simulators/utils/UnsupportedFlowKeywords.cpp index c8f819c03..58a144d11 100644 --- a/opm/simulators/utils/UnsupportedFlowKeywords.cpp +++ b/opm/simulators/utils/UnsupportedFlowKeywords.cpp @@ -398,7 +398,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords() {"NARROW", {false, std::nullopt}}, {"NCONSUMP", {false, std::nullopt}}, {"NEFAC", {false, std::nullopt}}, - {"NETBALAN", {false, std::nullopt}}, {"NETCOMPA", {false, std::nullopt}}, {"NEXT", {false, std::nullopt}}, {"NEXTSTPL", {false, std::nullopt}}, diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 15bc0ab7d..88c842474 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -282,7 +282,8 @@ namespace Opm { // at the beginning of each time step (Not report step) void prepareTimeStep(DeferredLogger& deferred_logger); void initPrimaryVariablesEvaluation() const; - void updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls); + bool shouldBalanceNetwork(const int reportStepIndex, const int iterationIdx) const; + std::pair updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls); void updateAndCommunicate(const int reportStepIdx, const int iterationIdx, @@ -369,6 +370,10 @@ namespace Opm { // and in the well equations. void assemble(const int iterationIdx, const double dt); + void assembleImpl(const int iterationIdx, + const double dt, + const std::size_t recursion_level, + DeferredLogger& local_deferredLogger); // called at the end of a time step void timeStepSucceeded(const double& simulationTime, const double dt); diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index b9cfb3913..266b51d11 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -2102,14 +2102,14 @@ inferLocalShutWells() } } -void +std::pair BlackoilWellModelGeneric:: updateNetworkPressures(const int reportStepIdx) { // Get the network and return if inactive. const auto& network = schedule()[reportStepIdx].network(); if (!network.active()) { - return; + return { false, 0.0 }; } node_pressures_ = WellGroupHelpers::computeNetworkPressures(network, this->wellState(), @@ -2119,6 +2119,8 @@ updateNetworkPressures(const int reportStepIdx) reportStepIdx); // Set the thp limits of wells + bool active_limit_change = false; + double network_imbalance = 0.0; for (auto& well : well_container_generic_) { // Producers only, since we so far only support the // "extended" network model (properties defined by @@ -2128,10 +2130,19 @@ updateNetworkPressures(const int reportStepIdx) 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. - well->setDynamicThpLimit(it->second); + const double new_limit = it->second; + well->setDynamicThpLimit(new_limit); + const SingleWellState& ws = this->wellState()[well->indexOfWell()]; + const bool thp_is_limit = ws.production_cmode == Well::ProducerCMode::THP; + const bool will_switch_to_thp = ws.thp < new_limit; + if (thp_is_limit || will_switch_to_thp) { + active_limit_change = true; + network_imbalance = std::max(network_imbalance, std::fabs(new_limit - ws.thp)); + } } } } + return { active_limit_change, network_imbalance }; } void diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index 56433b863..e7cc7e0cf 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -250,7 +250,7 @@ protected: bool wasDynamicallyShutThisTimeStep(const int well_index) const; - void updateNetworkPressures(const int reportStepIdx); + std::pair updateNetworkPressures(const int reportStepIdx); void updateWsolvent(const Group& group, const int reportStepIdx, diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 5034b28eb..5c8ad9972 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -870,7 +870,6 @@ namespace Opm { return; } - updatePerforationIntensiveQuantities(); if (iterationIdx == 0) { @@ -886,7 +885,27 @@ namespace Opm { OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed (It=0): ", terminal_output_, grid().comm()); } - updateWellControls(local_deferredLogger, /* check group controls */ true); + + assembleImpl(iterationIdx, dt, 0, local_deferredLogger); + + last_report_.converged = true; + last_report_.assemble_time_well += perfTimer.stop(); + } + + + + + + template + void + BlackoilWellModel:: + assembleImpl(const int iterationIdx, + const double dt, + const std::size_t recursion_level, + DeferredLogger& local_deferredLogger) + { + + const auto [network_changed, network_imbalance] = updateWellControls(local_deferredLogger, /* check group controls */ true); bool alq_updated = false; OPM_BEGIN_PARALLEL_TRY_CATCH(); @@ -911,10 +930,23 @@ namespace Opm { WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime, this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger); } - last_report_.converged = true; - last_report_.assemble_time_well += perfTimer.stop(); + + // Maybe do a recursive call to iterate network and well controls. + if (network_changed) { + if (shouldBalanceNetwork(reportStepIdx, iterationIdx)) { + const auto& balance = schedule()[reportStepIdx].network_balance(); + // Iterate if not converged, and number of iterations is not yet max (NETBALAN item 3). + if (recursion_level < balance.pressure_max_iter() && network_imbalance > balance.pressure_tolerance()) { + assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger); + } + } + } } + + + + template bool BlackoilWellModel:: @@ -1433,21 +1465,49 @@ namespace Opm { template - void + bool + BlackoilWellModel:: + shouldBalanceNetwork(const int reportStepIdx, const int iterationIdx) const + { + 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::pair BlackoilWellModel:: updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls) { // Even if there are no wells active locally, we cannot // return as the DeferredLogger uses global communication. // For no well active globally we simply return. - if( !wellsActive() ) return ; + if( !wellsActive() ) return { false, 0.0 }; const int episodeIdx = ebosSimulator_.episodeIndex(); const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations(); const auto& comm = ebosSimulator_.vanguard().grid().comm(); updateAndCommunicateGroupData(episodeIdx, iterationIdx); - updateNetworkPressures(episodeIdx); + const auto [local_network_changed, local_network_imbalance] + = shouldBalanceNetwork(episodeIdx, iterationIdx) ? + updateNetworkPressures(episodeIdx) : std::make_pair(false, 0.0); + const bool network_changed = comm.sum(local_network_changed); + const double network_imbalance = comm.max(local_network_imbalance); std::set switched_wells; @@ -1500,6 +1560,8 @@ namespace Opm { // update wsolvent fraction for REIN wells const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx); updateWsolvent(fieldGroup, episodeIdx, this->nupcolWellState()); + + return { network_changed, network_imbalance }; } diff --git a/opm/simulators/wells/WellInterfaceGeneric.cpp b/opm/simulators/wells/WellInterfaceGeneric.cpp index 48aca1de5..5530bb2e2 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.cpp +++ b/opm/simulators/wells/WellInterfaceGeneric.cpp @@ -347,6 +347,11 @@ void WellInterfaceGeneric::setDynamicThpLimit(const double thp_limit) dynamic_thp_limit_ = thp_limit; } +std::optional WellInterfaceGeneric::getDynamicThpLimit() const +{ + return dynamic_thp_limit_; +} + void WellInterfaceGeneric::updatePerforatedCell(std::vector& is_cell_perforated) { diff --git a/opm/simulators/wells/WellInterfaceGeneric.hpp b/opm/simulators/wells/WellInterfaceGeneric.hpp index 651f5ef7c..368693ad6 100644 --- a/opm/simulators/wells/WellInterfaceGeneric.hpp +++ b/opm/simulators/wells/WellInterfaceGeneric.hpp @@ -101,6 +101,7 @@ public: void setRepRadiusPerfLength(); void setWsolvent(const double wsolvent); void setDynamicThpLimit(const double thp_limit); + std::optional getDynamicThpLimit() const; void updatePerforatedCell(std::vector& is_cell_perforated); /// Returns true if the well has one or more THP limits/constraints.