diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 0603cd7263..b09d28efde 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -243,7 +243,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/wells/WellFilterCake.cpp opm/simulators/wells/WellGroupConstraints.cpp opm/simulators/wells/WellGroupControls.cpp - opm/simulators/wells/WellGroupHelpers.cpp + opm/simulators/wells/WellGroupHelper.cpp opm/simulators/wells/WellHelpers.cpp opm/simulators/wells/WellInterfaceFluidSystem.cpp opm/simulators/wells/WellInterfaceGeneric.cpp @@ -1159,7 +1159,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/WellFilterCake.hpp opm/simulators/wells/WellGroupConstraints.hpp opm/simulators/wells/WellGroupControls.hpp - opm/simulators/wells/WellGroupHelpers.hpp + opm/simulators/wells/WellGroupHelper.hpp opm/simulators/wells/WellHelpers.hpp opm/simulators/wells/WellInterfaceFluidSystem.hpp opm/simulators/wells/WellInterfaceGeneric.hpp diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 9348925795..9e9492fa8e 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -72,7 +72,7 @@ #include #include #include -#include +#include #include #include #include @@ -115,8 +115,8 @@ template class WellContributions; using ModelParameters = BlackoilModelParameters; using WellConnectionModule = WellConnectionAuxiliaryModule>; - using IndexTraits = typename FluidSystem::IndexTraitsType; + using WellGroupHelperType = WellGroupHelper; constexpr static std::size_t pressureVarIndex = GetPropType::pressureSwitchIdx; @@ -143,8 +143,6 @@ template class WellContributions; using AverageRegionalPressureType = RegionAverageCalculator:: AverageRegionalPressure >; - using WellGroupHelpersType = WellGroupHelpers; - explicit BlackoilWellModel(Simulator& simulator); void init(); diff --git a/opm/simulators/wells/BlackoilWellModelConstraints.cpp b/opm/simulators/wells/BlackoilWellModelConstraints.cpp index 4ec89cc128..25fb60de4c 100644 --- a/opm/simulators/wells/BlackoilWellModelConstraints.cpp +++ b/opm/simulators/wells/BlackoilWellModelConstraints.cpp @@ -30,7 +30,6 @@ #include #include -#include #include #include @@ -46,7 +45,6 @@ checkGroupInjectionConstraints(const Group& group, const int reportStepIdx, const Phase& phase) const { - const auto& well_state = wellModel_.wellState(); const auto& pu = wellModel_.phaseUsage(); int phasePos; @@ -65,13 +63,7 @@ checkGroupInjectionConstraints(const Group& group, if (currentControl != Group::InjectionCMode::RATE) { Scalar current_rate = 0.0; - current_rate += WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - phasePos, - /*isInjector*/true, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellSurfaceRates(group, phasePos, /*isInjector*/true); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -95,13 +87,7 @@ checkGroupInjectionConstraints(const Group& group, if (currentControl != Group::InjectionCMode::RESV) { Scalar current_rate = 0.0; - current_rate += WellGroupHelpersType::sumWellResRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - phasePos, - /*isInjector*/true, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellResRates(group, phasePos, /*is_injector=*/true); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -126,25 +112,13 @@ checkGroupInjectionConstraints(const Group& group, Scalar production_Rate = 0.0; const auto& controls = group.injectionControls(phase, wellModel_.summaryState()); const Group& groupRein = wellModel_.schedule().getGroup(controls.reinj_group, reportStepIdx); - production_Rate += WellGroupHelpersType::sumWellSurfaceRates(groupRein, - wellModel_.schedule(), - well_state, - reportStepIdx, - phasePos, - /*isInjector*/false, - wellModel_.summaryState()); + production_Rate += wgHelper().sumWellSurfaceRates(groupRein, phasePos, /*is_injector=*/false); // sum over all nodes production_Rate = wellModel_.comm().sum(production_Rate); Scalar current_rate = 0.0; - current_rate += WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - phasePos, - /*isInjector*/true, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellSurfaceRates(group, phasePos, /*is_injector=*/true); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -164,52 +138,29 @@ checkGroupInjectionConstraints(const Group& group, Scalar voidage_rate = 0.0; const auto& controls = group.injectionControls(phase, wellModel_.summaryState()); const Group& groupVoidage = wellModel_.schedule().getGroup(controls.voidage_group, reportStepIdx); - voidage_rate += WellGroupHelpersType::sumWellResRates(groupVoidage, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(waterPhaseIdx), - false, - wellModel_.summaryState()); - voidage_rate += WellGroupHelpersType::sumWellResRates(groupVoidage, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(oilPhaseIdx), - false, - wellModel_.summaryState()); - voidage_rate += WellGroupHelpersType::sumWellResRates(groupVoidage, - wellModel_.schedule(), - well_state, reportStepIdx, - pu.canonicalToActivePhaseIdx(gasPhaseIdx), - false, - wellModel_.summaryState()); + voidage_rate += wgHelper().sumWellResRates(groupVoidage, + pu.canonicalToActivePhaseIdx(waterPhaseIdx), + /*is_injector=*/false); + voidage_rate += wgHelper().sumWellResRates(groupVoidage, + pu.canonicalToActivePhaseIdx(oilPhaseIdx), + /*is_injector=*/false); + voidage_rate += wgHelper().sumWellResRates(groupVoidage, + pu.canonicalToActivePhaseIdx(gasPhaseIdx), + /*is_injector=*/false); // sum over all nodes voidage_rate = wellModel_.comm().sum(voidage_rate); Scalar total_rate = 0.0; - total_rate += WellGroupHelpersType::sumWellResRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(waterPhaseIdx), - true, - wellModel_.summaryState()); - total_rate += WellGroupHelpersType::sumWellResRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(oilPhaseIdx), - true, - wellModel_.summaryState()); - total_rate += WellGroupHelpersType::sumWellResRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(gasPhaseIdx), - true, - wellModel_.summaryState()); + total_rate += wgHelper().sumWellResRates(group, + pu.canonicalToActivePhaseIdx(waterPhaseIdx), + /*is_injector=*/true); + total_rate += wgHelper().sumWellResRates(group, + pu.canonicalToActivePhaseIdx(oilPhaseIdx), + /*is_injector=*/true); + total_rate += wgHelper().sumWellResRates(group, + pu.canonicalToActivePhaseIdx(gasPhaseIdx), + /*is_injector=*/true); // sum over all nodes total_rate = wellModel_.comm().sum(total_rate); @@ -228,12 +179,8 @@ checkGroupInjectionConstraints(const Group& group, template std::pair BlackoilWellModelConstraints:: -checkGroupProductionConstraints(const Group& group, - const int reportStepIdx, - DeferredLogger& deferred_logger) const +checkGroupProductionConstraints(const Group& group, DeferredLogger& deferred_logger) const { - const auto& well_state = wellModel_.wellState(); - const auto controls = group.productionControls(wellModel_.summaryState()); const Group::ProductionCMode& currentControl = wellModel_.groupState().production_control(group.name()); const auto& pu = wellModel_.phaseUsage(); @@ -242,13 +189,9 @@ checkGroupProductionConstraints(const Group& group, if (currentControl != Group::ProductionCMode::ORAT) { Scalar current_rate = 0.0; - current_rate += WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(oilPhaseIdx), - false, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellSurfaceRates(group, + pu.canonicalToActivePhaseIdx(oilPhaseIdx), + /*is_injector=*/false); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -267,13 +210,9 @@ checkGroupProductionConstraints(const Group& group, if (currentControl != Group::ProductionCMode::WRAT) { Scalar current_rate = 0.0; - current_rate += WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(waterPhaseIdx), - false, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellSurfaceRates(group, + pu.canonicalToActivePhaseIdx(waterPhaseIdx), + /*is_injector=*/false); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -291,13 +230,9 @@ checkGroupProductionConstraints(const Group& group, if (currentControl != Group::ProductionCMode::GRAT) { Scalar current_rate = 0.0; - current_rate += WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(gasPhaseIdx), - false, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellSurfaceRates(group, + pu.canonicalToActivePhaseIdx(gasPhaseIdx), + /*is_injector=*/false); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -314,20 +249,12 @@ checkGroupProductionConstraints(const Group& group, if (currentControl != Group::ProductionCMode::LRAT) { Scalar current_rate = 0.0; - current_rate += WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(oilPhaseIdx), - false, - wellModel_.summaryState()); - current_rate += WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(waterPhaseIdx), - false, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellSurfaceRates(group, + pu.canonicalToActivePhaseIdx(oilPhaseIdx), + /*is_injector=*/false); + current_rate += wgHelper().sumWellSurfaceRates(group, + pu.canonicalToActivePhaseIdx(waterPhaseIdx), + /*is_injector=*/false); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -335,13 +262,9 @@ checkGroupProductionConstraints(const Group& group, bool skip = false; if (controls.liquid_target == controls.oil_target) { Scalar current_water_rate = - WellGroupHelpersType::sumWellSurfaceRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(waterPhaseIdx), - false, - wellModel_.summaryState()); + wgHelper().sumWellSurfaceRates(group, + pu.canonicalToActivePhaseIdx(waterPhaseIdx), + /*is_injector=*/false); current_water_rate = wellModel_.comm().sum(current_water_rate); if (std::abs(current_water_rate) < 1e-12) { skip = true; @@ -367,30 +290,15 @@ checkGroupProductionConstraints(const Group& group, if (currentControl != Group::ProductionCMode::RESV) { Scalar current_rate = 0.0; - current_rate += - WellGroupHelpersType::sumWellResRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(waterPhaseIdx), - false, - wellModel_.summaryState()); - current_rate += - WellGroupHelpersType::sumWellResRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(oilPhaseIdx), - false, - wellModel_.summaryState()); - current_rate += - WellGroupHelpersType::sumWellResRates(group, - wellModel_.schedule(), - well_state, - reportStepIdx, - pu.canonicalToActivePhaseIdx(gasPhaseIdx), - false, - wellModel_.summaryState()); + current_rate += wgHelper().sumWellResRates(group, + pu.canonicalToActivePhaseIdx(waterPhaseIdx), + /*is_injector=*/false); + current_rate += wgHelper().sumWellResRates(group, + pu.canonicalToActivePhaseIdx(oilPhaseIdx), + /*is_injector=*/false); + current_rate += wgHelper().sumWellResRates(group, + pu.canonicalToActivePhaseIdx(gasPhaseIdx), + /*is_injector=*/false); // sum over all nodes current_rate = wellModel_.comm().sum(current_rate); @@ -434,9 +342,7 @@ checkGroupConstraints(const Group& group, } } if (group.isProductionGroup()) { - const auto& check = this->checkGroupProductionConstraints(group, - reportStepIdx, - deferred_logger); + const auto& check = this->checkGroupProductionConstraints(group, deferred_logger); if (check.first != Group::ProductionCMode::NONE) { return true; @@ -478,10 +384,8 @@ actionOnBrokenConstraints(const Group& group, template bool BlackoilWellModelConstraints:: actionOnBrokenConstraints(const Group& group, - const int reportStepIdx, const Group::GroupLimitAction group_limit_action, const Group::ProductionCMode& newControl, - const WellState& well_state, std::optional& worst_offending_well, GroupState& group_state, DeferredLogger& deferred_logger) const @@ -548,9 +452,7 @@ actionOnBrokenConstraints(const Group& group, } case Group::ExceedAction::WELL: { std::tie(worst_offending_well, std::ignore) = - WellGroupHelpersType::worstOffendingWell(group, wellModel_.schedule(), reportStepIdx, - newControl, - wellModel_.comm(), well_state, deferred_logger); + wgHelper().worstOffendingWell(group, newControl, wellModel_.comm(), deferred_logger); break; } case Group::ExceedAction::PLUG: { @@ -642,13 +544,10 @@ updateGroupIndividualControl(const Group& group, } this->actionOnBrokenConstraints(group, changed_this.first, phase, group_state, deferred_logger); - WellGroupHelpersType::updateWellRatesFromGroupTargetScale(changed_this.second, - group, - wellModel_.schedule(), - reportStepIdx, - /* isInjector */ false, - wellModel_.groupState(), - well_state); + wgHelper().updateWellRatesFromGroupTargetScale(changed_this.second, + group, + /*is_injector=*/false, + well_state); changed = true; } } @@ -676,32 +575,27 @@ updateGroupIndividualControl(const Group& group, } } - const auto& changed_this = this->checkGroupProductionConstraints(group, - reportStepIdx, - deferred_logger); + const auto& changed_this = this->checkGroupProductionConstraints(group, deferred_logger); const auto controls = group.productionControls(wellModel_.summaryState()); if (changed_this.first != Group::ProductionCMode::NONE) { std::optional worst_offending_well = std::nullopt; - changed = this->actionOnBrokenConstraints(group, reportStepIdx, - controls.group_limit_action, - changed_this.first, well_state, - worst_offending_well, - group_state, deferred_logger); + changed = this->actionOnBrokenConstraints(group, + controls.group_limit_action, + changed_this.first, + worst_offending_well, + group_state, deferred_logger); if(changed) { if (update_group_switching_log || switched_prod[group.name()].empty()) { switched_prod[group.name()].push_back(currentControl); } - WellGroupHelpersType::updateWellRatesFromGroupTargetScale(changed_this.second, - group, - wellModel_.schedule(), - reportStepIdx, - /* isInjector */ false, - wellModel_.groupState(), - well_state); + wgHelper().updateWellRatesFromGroupTargetScale(changed_this.second, + group, + /*is_injector=*/false, + well_state); } else if (worst_offending_well) { closed_offending_wells.insert_or_assign(group.name(), std::make_pair(Group::ProductionCMode2String(changed_this.first), *worst_offending_well)); diff --git a/opm/simulators/wells/BlackoilWellModelConstraints.hpp b/opm/simulators/wells/BlackoilWellModelConstraints.hpp index 3ee5b1942f..d9b5a3a862 100644 --- a/opm/simulators/wells/BlackoilWellModelConstraints.hpp +++ b/opm/simulators/wells/BlackoilWellModelConstraints.hpp @@ -34,7 +34,7 @@ class DeferredLogger; template class GroupState; class SummaryState; template class WellState; -template class WellGroupHelpers; +template class WellGroupHelper; /// Class for handling constraints for the blackoil well model. template @@ -46,7 +46,7 @@ public: constexpr static int oilPhaseIdx = IndexTraits::oilPhaseIdx; constexpr static int gasPhaseIdx = IndexTraits::gasPhaseIdx; - using WellGroupHelpersType = WellGroupHelpers; + using WellGroupHelperType = WellGroupHelper; //! \brief Constructor initializes reference to the well model. explicit BlackoilWellModelConstraints(const BlackoilWellModelGeneric& wellModel) @@ -67,10 +67,8 @@ public: //! \brief Execute action on broken constraint for a production well group. Return true if a group control is changed bool actionOnBrokenConstraints(const Group& group, - const int reportStepIdx, const Group::GroupLimitAction group_limit_action, const Group::ProductionCMode& newControl, - const WellState& well_state, std::optional& worst_offending_well, GroupState& group_state, DeferredLogger& deferred_logger) const; @@ -97,9 +95,9 @@ private: //! \brief Check and return value and type of constraints for a production well group. std::pair checkGroupProductionConstraints(const Group& group, - const int reportStepIdx, DeferredLogger& deferred_logger) const; + const WellGroupHelperType& wgHelper() const { return wellModel_.wgHelper(); } const BlackoilWellModelGeneric& wellModel_; //!< Reference to well model }; diff --git a/opm/simulators/wells/BlackoilWellModelGasLift_impl.hpp b/opm/simulators/wells/BlackoilWellModelGasLift_impl.hpp index 7b05b9cd6b..ff6e5aa773 100644 --- a/opm/simulators/wells/BlackoilWellModelGasLift_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModelGasLift_impl.hpp @@ -334,7 +334,8 @@ updateWellPotentials(const Simulator& simulator, auto cur_exc_type = ExceptionType::NONE; try { std::vector potentials; - well->computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger); + const auto& wgHelper = simulator.problem().wellModel().wgHelper(); + well->computeWellPotentials(simulator, well_state_copy, wgHelper, potentials, deferred_logger); auto& ws = wellState.well(well->indexOfWell()); for (int p = 0; p < np; ++p) { // make sure the potentials are positive diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 2fe3e8cf59..fef06c75f1 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -64,7 +64,7 @@ #include #include #include -#include +#include #include #include @@ -109,6 +109,7 @@ BlackoilWellModelGeneric(Schedule& schedule, , active_wgstate_(pu) , last_valid_wgstate_(pu) , nupcol_wgstate_(pu) + , wg_helper_(this->wellState(), this->groupState(), this->schedule(), summaryState, guideRate_, pu) { const auto numProcs = comm_.size(); @@ -505,20 +506,12 @@ checkGconsaleLimits(const Group& group, const Group::ProductionCMode& oldProductionControl = this->groupState().production_control(group.name()); const int gasPos = phase_usage_info_.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); - Scalar production_rate = WellGroupHelpersType::sumWellSurfaceRates(group, - schedule(), - well_state, - reportStepIdx, - gasPos, - /*isInjector*/false, - this->summaryState_); - Scalar injection_rate = WellGroupHelpersType::sumWellSurfaceRates(group, - schedule(), - well_state, - reportStepIdx, - gasPos, - /*isInjector*/true, - this->summaryState_); + Scalar production_rate = this->wgHelper().sumWellSurfaceRates( + group, gasPos, /*is_injector=*/false + ); + Scalar injection_rate = this->wgHelper().sumWellSurfaceRates( + group, gasPos, /*is_injector=*/true + ); // sum over all nodes injection_rate = comm_.sum(injection_rate); production_rate = comm_.sum(production_rate); @@ -655,13 +648,9 @@ checkGroupHigherConstraints(const Group& group, const std::vector reduction_rates = this->groupState().injection_reduction_rates(group.name()); for (int phasePos = 0; phasePos < this->numPhases(); ++phasePos) { - const Scalar local_current_rate = WellGroupHelpersType::sumWellSurfaceRates(group, - schedule(), - this->wellState(), - reportStepIdx, - phasePos, - /* isInjector */ true, - this->summaryState_); + const Scalar local_current_rate = this->wgHelper().sumWellSurfaceRates( + group, phasePos, /*is_injector=*/true + ); // Sum over all processes rates[phasePos] = comm_.sum(local_current_rate) - reduction_rates[phasePos]; } @@ -698,22 +687,17 @@ checkGroupHigherConstraints(const Group& group, // Check higher up only if under individual (not FLD) control. if (currentControl != Group::InjectionCMode::FLD && group.injectionGroupControlAvailable(phase)) { const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx); - const auto [is_changed, scaling_factor] = - WellGroupHelpersType::checkGroupConstraintsInj(group.name(), - group.parent(), - parentGroup, - this->wellState(), - this->groupState(), - reportStepIdx, - &guideRate_, - rates.data(), - phase, - group.getGroupEfficiencyFactor(), - schedule(), - summaryState_, - resv_coeff_inj, - /*check_guide_rate*/true, - deferred_logger); + const auto [is_changed, scaling_factor] = this->wgHelper().checkGroupConstraintsInj( + group.name(), + group.parent(), + parentGroup, + rates.data(), + phase, + group.getGroupEfficiencyFactor(), + resv_coeff_inj, + /*check_guide_rate*/true, + deferred_logger + ); if (is_changed) { auto& group_log = switched_inj_groups_[group.name()][static_cast>(phase)]; if (update_group_switching_log || group_log.empty()) { @@ -723,13 +707,9 @@ checkGroupHigherConstraints(const Group& group, actionOnBrokenConstraints(group, Group::InjectionCMode::FLD, phase, this->groupState(), deferred_logger); - WellGroupHelpersType::updateWellRatesFromGroupTargetScale(scaling_factor, - group, - schedule(), - reportStepIdx, - /* isInjector */ true, - this->groupState(), - this->wellState()); + this->wgHelper().updateWellRatesFromGroupTargetScale( + scaling_factor, group, /*is_injector=*/true, this->wellState() + ); changed = true; } } @@ -762,13 +742,9 @@ checkGroupHigherConstraints(const Group& group, } } for (int phasePos = 0; phasePos < this->numPhases(); ++phasePos) { - const Scalar local_current_rate = WellGroupHelpersType::sumWellSurfaceRates(group, - schedule(), - this->wellState(), - reportStepIdx, - phasePos, - /* isInjector */ false, - this->summaryState_); + const Scalar local_current_rate = this->wgHelper().sumWellSurfaceRates( + group, phasePos, /*is_injector=*/false + ); // Sum over all processes rates[phasePos] = -comm_.sum(local_current_rate) - reduction_rates[phasePos]; } @@ -777,28 +753,23 @@ checkGroupHigherConstraints(const Group& group, // Check higher up only if under individual (not FLD) control. if (currentControl != Group::ProductionCMode::FLD && group.productionGroupControlAvailable()) { const Group& parentGroup = schedule().getGroup(group.parent(), reportStepIdx); - const auto [is_changed, scaling_factor] = - WellGroupHelpersType::checkGroupConstraintsProd(group.name(), - group.parent(), - parentGroup, - this->wellState(), - this->groupState(), - reportStepIdx, - &guideRate_, - rates.data(), - group.getGroupEfficiencyFactor(), - schedule(), - summaryState_, - resv_coeff, - /*check_guide_rate*/true, - deferred_logger); + const auto [is_changed, scaling_factor] = this->wgHelper().checkGroupConstraintsProd( + group.name(), + group.parent(), + parentGroup, + rates.data(), + group.getGroupEfficiencyFactor(), + resv_coeff, + /*check_guide_rate*/true, + deferred_logger + ); if (is_changed) { const auto group_limit_action = group.productionControls(summaryState_).group_limit_action; std::optional worst_offending_well = std::nullopt; changed = BlackoilWellModelConstraints(*this). - actionOnBrokenConstraints(group, reportStepIdx, group_limit_action, + actionOnBrokenConstraints(group, + group_limit_action, Group::ProductionCMode::FLD, - this->wellState(), worst_offending_well, this->groupState(), deferred_logger); @@ -807,13 +778,9 @@ checkGroupHigherConstraints(const Group& group, if (update_group_switching_log || switched_prod_groups_[group.name()].empty()) { switched_prod_groups_[group.name()].push_back(currentControl); } - WellGroupHelpersType::updateWellRatesFromGroupTargetScale(scaling_factor, - group, - schedule(), - reportStepIdx, - /* isInjector */ false, - this->groupState(), - this->wellState()); + this->wgHelper().updateWellRatesFromGroupTargetScale( + scaling_factor, group, /*is_injector=*/false, this->wellState() + ); } } } @@ -1006,18 +973,12 @@ updateWsolvent(const Group& group, const int gasPos = phase_usage_info_.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); const auto& controls = group.injectionControls(Phase::GAS, summaryState_); const Group& groupRein = schedule_.getGroup(controls.reinj_group, reportStepIdx); - Scalar gasProductionRate = WellGroupHelpersType::sumWellSurfaceRates(groupRein, - schedule_, - wellState, - reportStepIdx, - gasPos, - /*isInjector*/false, - this->summaryState_); - Scalar solventProductionRate = WellGroupHelpersType::sumSolventRates(groupRein, - schedule_, - wellState, - reportStepIdx, - /*isInjector*/false); + Scalar gasProductionRate = this->wgHelper().sumWellSurfaceRates( + groupRein, gasPos, /*is_injector=*/false + ); + Scalar solventProductionRate = this->wgHelper().sumSolventRates( + groupRein, /*is_injector=*/false + ); solventProductionRate = comm_.sum(solventProductionRate); gasProductionRate = comm_.sum(gasProductionRate); @@ -1260,13 +1221,11 @@ assignNodeValues(std::map& nodevalues, return; } - auto converged_pressures = WellGroupHelpersType::computeNetworkPressures(network, - this->wellState(), - this->groupState(), - *(vfp_properties_->getProd()), - schedule(), - comm_, - reportStepIdx); + auto converged_pressures = this->wgHelper().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() ); @@ -1327,6 +1286,7 @@ updateAndCommunicateGroupData(const int reportStepIdx, this->wellState().updateGlobalIsGrup(comm_, well_status); + WellGroupHelperType &wg_helper = this->wgHelper(); if (iterationIdx < nupcol) { OPM_TIMEBLOCK(updateNupcol); this->updateNupcolWGState(); @@ -1342,26 +1302,27 @@ updateAndCommunicateGroupData(const int reportStepIdx, const Group& group = schedule().getGroup(gr_name, reportStepIdx); const int np = this->wellState().numPhases(); Scalar gr_rate_nupcol = 0.0; - for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - gr_rate_nupcol += WellGroupHelpersType::sumWellPhaseRates(is_vrep, - group, - schedule(), - this->nupcolWellState(), - this->summaryState_, - reportStepIdx, - phaseIdx, - /*isInjector*/ false); + { + // Temporarily use the nupcol well state for all helper functions + // At the end of this scope, the well state will be restored to its original value + auto guard = wg_helper.pushWellState(this->nupcolWellState()); + for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { + gr_rate_nupcol += wg_helper.sumWellPhaseRates( + /*res_rates=*/is_vrep, + group, + phaseIdx, + /*is_injector=*/false + ); + } } Scalar gr_rate = 0.0; for (int phaseIdx = 0; phaseIdx < np; ++phaseIdx) { - gr_rate += WellGroupHelpersType::sumWellPhaseRates(is_vrep, - group, - schedule(), - this->wellState(), - this->summaryState_, - reportStepIdx, - phaseIdx, - /*isInjector*/ false); + gr_rate += wg_helper.sumWellPhaseRates( + /*res_rates=*/is_vrep, + group, + phaseIdx, + /*is_injector=*/false + ); } // sum contributions from owned wells to everybody gr_rate_nupcol = comm_.sum(gr_rate_nupcol); @@ -1385,93 +1346,33 @@ updateAndCommunicateGroupData(const int reportStepIdx, } } } - auto& well_state = this->wellState(); - const auto& well_state_nupcol = this->nupcolWellState(); - - constexpr int num_configs = 4; - constexpr std::array is_production_group = {true, false, false, false}; - constexpr std::array phases = { Phase::OIL, Phase::WATER, Phase::OIL, Phase::GAS }; - for (int i = 0; i < num_configs; i++) { - WellGroupHelpersType::updateGroupControlledWells(schedule(), - well_state, - this->groupState(), - summaryState_, - &guideRate_, - reportStepIdx, - "FIELD", - is_production_group[i], - phases[i]); + { + constexpr int num_configs = 4; + constexpr std::array is_production_group = {true, false, false, false}; + constexpr std::array phases = { Phase::OIL, Phase::WATER, Phase::OIL, Phase::GAS }; + for (int i = 0; i < num_configs; i++) { + wg_helper.updateGroupControlledWells( + is_production_group[i], phases[i], this->groupState(), deferred_logger + ); + } } - // the group target reduction rates needs to be update since wells may have switched to/from GRUP control // The group target reduction does not honor NUPCOL. - std::vector groupTargetReduction(numPhases(), 0.0); - WellGroupHelpersType::updateGroupTargetReduction(fieldGroup, - schedule(), - reportStepIdx, - /*isInjector*/ false, - guideRate_, - well_state, - summaryState_, - this->groupState(), - groupTargetReduction); - std::vector groupTargetReductionInj(numPhases(), 0.0); - WellGroupHelpersType::updateGroupTargetReduction(fieldGroup, - schedule(), - reportStepIdx, - /*isInjector*/ true, - guideRate_, - well_state, - summaryState_, - this->groupState(), - groupTargetReductionInj); - - WellGroupHelpersType::updateREINForGroups(fieldGroup, - schedule(), - reportStepIdx, - summaryState_, - well_state_nupcol, - this->groupState(), - comm_.rank() == 0); - WellGroupHelpersType::updateVREPForGroups(fieldGroup, - schedule(), - reportStepIdx, - well_state_nupcol, - this->groupState(), - this->summaryState_); - - WellGroupHelpersType::updateReservoirRatesInjectionGroups(fieldGroup, - schedule(), - reportStepIdx, - well_state_nupcol, - this->groupState(), - this->summaryState_); - WellGroupHelpersType::updateSurfaceRatesInjectionGroups(fieldGroup, - schedule(), - reportStepIdx, - well_state_nupcol, - this->groupState(), - this->summaryState_); - WellGroupHelpersType::updateNetworkLeafNodeProductionRates(schedule(), - reportStepIdx, - well_state_nupcol, - this->groupState(), - this->summaryState_); - - WellGroupHelpersType::updateGroupProductionRates(fieldGroup, - schedule(), - reportStepIdx, - well_state_nupcol, - this->groupState(), - this->summaryState_); - - WellGroupHelpersType::updateWellRates(fieldGroup, - schedule(), - reportStepIdx, - well_state_nupcol, - well_state); - - well_state.communicateGroupRates(comm_); + wg_helper.updateGroupTargetReduction(fieldGroup, /*is_injector=*/false, this->groupState()); + wg_helper.updateGroupTargetReduction(fieldGroup, /*is_injector=*/true, this->groupState()); + { + // Temporarily use the nupcol well state for all helper functions + // At the end of this scope, the well state will be restored to its original value + auto guard = wg_helper.pushWellState(this->nupcolWellState()); + wg_helper.updateREINForGroups(fieldGroup, /*sum_rank=*/comm_.rank() == 0, this->groupState()); + wg_helper.updateVREPForGroups(fieldGroup, this->groupState()); + wg_helper.updateReservoirRatesInjectionGroups(fieldGroup, this->groupState()); + wg_helper.updateSurfaceRatesInjectionGroups(fieldGroup, this->groupState()); + wg_helper.updateNetworkLeafNodeProductionRates(this->groupState()); + wg_helper.updateGroupProductionRates(fieldGroup, this->groupState()); + } + wg_helper.updateWellRates(fieldGroup, this->nupcolWellState(), this->wellState()); + this->wellState().communicateGroupRates(comm_); this->groupState().communicate_rates(comm_); if (update_wellgrouptarget) { @@ -1487,20 +1388,15 @@ updateAndCommunicateGroupData(const int reportStepIdx, // Translate injector type from control to Phase. Scalar group_target = std::numeric_limits::max(); if (well->isProducer()) { - group_target = - WellGroupHelpersType::getWellGroupTargetProducer(well->name(), - well->wellEcl().groupName(), - group, - this->wellState(), - this->groupState(), - well->currentStep(), - well->guideRate(), - ws.surface_rates.data(), - efficiencyFactor, - this->schedule(), - summaryState_, - resv_coeff, - deferred_logger); + group_target = wg_helper.getWellGroupTargetProducer( + well->name(), + well->wellEcl().groupName(), + group, + ws.surface_rates.data(), + efficiencyFactor, + resv_coeff, + deferred_logger + ); } else { const auto& well_controls = well->wellEcl().injectionControls(summaryState_); auto injectorType = well_controls.injector_type; @@ -1524,21 +1420,16 @@ updateAndCommunicateGroupData(const int reportStepIdx, default: throw std::logic_error("MULTI-phase injection is not supported, but was requested for well " + well->name()); } - group_target = - WellGroupHelpersType::getWellGroupTargetInjector(well->name(), - well->wellEcl().groupName(), - group, - this->wellState(), - this->groupState(), - well->currentStep(), - well->guideRate(), - ws.surface_rates.data(), - injectionPhase, - efficiencyFactor, - this->schedule(), - summaryState_, - resv_coeff, - deferred_logger); + group_target = wg_helper.getWellGroupTargetInjector( + well->name(), + well->wellEcl().groupName(), + group, + ws.surface_rates.data(), + injectionPhase, + efficiencyFactor, + resv_coeff, + deferred_logger + ); } auto& ws_update = this->wellState().well(well->indexOfWell()); ws_update.group_target = group_target; @@ -1675,13 +1566,8 @@ updateNetworkPressures(const int reportStepIdx, const Scalar damping_factor, con const auto previous_node_pressures = node_pressures_; - node_pressures_ = WellGroupHelpersType::computeNetworkPressures(network, - this->wellState(), - this->groupState(), - *(vfp_properties_->getProd()), - schedule(), - comm_, - reportStepIdx); + node_pressures_ = this->wgHelper().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.; @@ -1748,11 +1634,10 @@ calculateEfficiencyFactors(const int reportStepIdx) const Well& wellEcl = well->wellEcl(); Scalar well_efficiency_factor = wellEcl.getEfficiencyFactor() * wellState().getGlobalEfficiencyScalingFactor(well->name()); - WellGroupHelpersType::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), - reportStepIdx), - schedule(), - reportStepIdx, - well_efficiency_factor); + this->wgHelper().accumulateGroupEfficiencyFactor( + schedule().getGroup(wellEcl.groupName(), reportStepIdx), + well_efficiency_factor + ); well->setWellEfficiencyFactor(well_efficiency_factor); } } diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp index f6c8ac31c9..9b546db66c 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp @@ -39,6 +39,7 @@ #include #include #include +#include #include #include #include @@ -69,7 +70,6 @@ namespace Opm { struct SimulatorUpdate; class SummaryConfig; template class VFPProperties; - template class WellGroupHelpers; template class WellInterfaceGeneric; template class WellState; } // namespace Opm @@ -93,7 +93,7 @@ namespace Opm { template class BlackoilWellModelGeneric { - using WellGroupHelpersType = WellGroupHelpers; + using WellGroupHelperType = WellGroupHelper; public: BlackoilWellModelGeneric(Schedule& schedule, BlackoilWellModelGasLiftGeneric& gaslift, @@ -165,11 +165,19 @@ public: /* Will return the currently active nupcolWellState; must update the internal nupcol wellstate with updateNupcolWGState() first. + + Both const and non-const accessors are provided. The non-const + accessor is required for the WellStateGuard pattern and pushWellState() + in WellGroupHelper, which temporarily switches WellGroupHelper to use this state. */ const WellState& nupcolWellState() const { return this->nupcol_wgstate_.well_state; } + WellState& nupcolWellState() + { + return this->nupcol_wgstate_.well_state; + } GroupState& groupState() { return this->active_wgstate_.group_state; } WellTestState& wellTestState() { return this->active_wgstate_.well_test_state; } @@ -295,6 +303,9 @@ public: const ConnectionIndexMap& connectionIndexMap(const std::size_t idx) { return conn_idx_map_[idx]; } + WellGroupHelperType& wgHelper() { return wg_helper_; } + const WellGroupHelperType& wgHelper() const { return wg_helper_; } + protected: /* The dynamic state of the well model is maintained with an instance @@ -352,6 +363,8 @@ protected: { 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->wg_helper_.updateState(this->wellState(), this->groupState()); } /* @@ -575,6 +588,7 @@ protected: WGState active_wgstate_; WGState last_valid_wgstate_; WGState nupcol_wgstate_; + WellGroupHelperType wg_helper_; WellGroupEvents report_step_start_events_; //!< Well group events at start of report step bool wellStructureChangedDynamically_{false}; diff --git a/opm/simulators/wells/BlackoilWellModelGuideRates.cpp b/opm/simulators/wells/BlackoilWellModelGuideRates.cpp index db63de4fd5..10f2d488ee 100644 --- a/opm/simulators/wells/BlackoilWellModelGuideRates.cpp +++ b/opm/simulators/wells/BlackoilWellModelGuideRates.cpp @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include @@ -331,8 +331,7 @@ getGuideRateValues(const Well& well) const return grval; } - const auto qs = WellGroupHelpers:: - getWellRateVector(wellModel_.wellState(), wname); + const auto qs = this->wgHelper().getWellRateVector(wname); this->getGuideRateValues(qs, well.isInjector(), wname, grval); @@ -359,8 +358,7 @@ getGuideRateValues(const Group& group) const return grval; } - const auto qs = WellGroupHelpers:: - getProductionGroupRateVector(wellModel_.groupState(), wellModel_.phaseUsage(), gname); + const auto qs = this->wgHelper().getProductionGroupRateVector(gname); const auto is_inj = false; // This procedure only applies to G*PGR. this->getGuideRateValues(qs, is_inj, gname, grval); @@ -439,8 +437,7 @@ calculateWellGuideRates(const int reportStepIdx) const const auto get_gr = parent || RetrieveWellGuideRate{wellModel_.guideRate(), wname}; - const auto qs = WellGroupHelpers:: - getWellRateVector(wellModel_.wellState(), wname); + const auto qs = this->wgHelper().getWellRateVector(wname); auto getGR = [this, &wname, &qs](const GuideRateModel::Target t) { diff --git a/opm/simulators/wells/BlackoilWellModelGuideRates.hpp b/opm/simulators/wells/BlackoilWellModelGuideRates.hpp index 86a7e8adec..6be3364137 100644 --- a/opm/simulators/wells/BlackoilWellModelGuideRates.hpp +++ b/opm/simulators/wells/BlackoilWellModelGuideRates.hpp @@ -30,6 +30,7 @@ namespace Opm { +template class WellGroupHelper; template class BlackoilWellModelGeneric; namespace data { struct GroupData; @@ -45,6 +46,8 @@ template class BlackoilWellModelGuideRates { public: + using WellGroupHelperType = WellGroupHelper; + //! \brief Constructor initializes reference to the well model. explicit BlackoilWellModelGuideRates(const BlackoilWellModelGeneric& wellModel) : wellModel_(wellModel) @@ -70,6 +73,8 @@ public: //! \brief Check if a guide rate update is needed. bool guideRateUpdateIsNeeded(const int reportStepIdx) const; + const WellGroupHelperType& wgHelper() const { return wellModel_.wgHelper(); } + private: //! \brief Obtain guide rate values. void getGuideRateValues(const GuideRate::RateVector& qs, diff --git a/opm/simulators/wells/BlackoilWellModelNldd_impl.hpp b/opm/simulators/wells/BlackoilWellModelNldd_impl.hpp index 2d04825d1c..d4f3ff71d0 100644 --- a/opm/simulators/wells/BlackoilWellModelNldd_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModelNldd_impl.hpp @@ -63,8 +63,8 @@ assembleWellEq(const double dt, if (this->well_domain().at(well->name()) == domain.index) { well->assembleWellEq(wellModel_.simulator(), dt, + wellModel_.wgHelper(), wellModel_.wellState(), - wellModel_.groupState(), deferred_logger); } } @@ -192,8 +192,8 @@ updateWellControls(DeferredLogger& deferred_logger, constexpr auto mode = WellInterface::IndividualOrGroup::Individual; well->updateWellControl(wellModel_.simulator(), mode, + wellModel_.wgHelper(), wellModel_.wellState(), - wellModel_.groupState(), deferred_logger); } } diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 0f97b65d28..2a095012fc 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -49,7 +49,7 @@ #include #include #include -#include +#include #include #include @@ -195,6 +195,7 @@ namespace Opm { { DeferredLogger local_deferredLogger{}; + this->wgHelper().setReportStep(timeStepIdx); this->report_step_starts_ = true; this->report_step_start_events_ = this->schedule()[timeStepIdx].wellgroup_events(); @@ -299,20 +300,16 @@ namespace Opm { const auto& fieldGroup = this->schedule().getGroup("FIELD", reportStepIdx); - WellGroupHelpersType::setCmodeGroup(fieldGroup, - this->schedule(), - this->summaryState(), - reportStepIdx, - this->groupState()); + this->wgHelper().setCmodeGroup(fieldGroup, this->groupState()); // Define per region average pressure calculators for use by // pressure maintenance groups (GPMAINT keyword). if (this->schedule()[reportStepIdx].has_gpmaint()) { - WellGroupHelpersType::setRegionAveragePressureCalculator(fieldGroup, - this->schedule(), - reportStepIdx, - this->eclState_.fieldProps(), - this->regionalAveragePressureCalculator_); + this->wgHelper().setRegionAveragePressureCalculator( + fieldGroup, + this->eclState_.fieldProps(), + this->regionalAveragePressureCalculator_ + ); } } OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, @@ -364,10 +361,8 @@ namespace Opm { } this->resetWGState(); - const int reportStepIdx = simulator_.episodeIndex(); - this->wellState().updateWellsDefaultALQ(this->schedule(), reportStepIdx, this->summaryState()); this->wellState().gliftTimeStepInit(); @@ -482,13 +477,10 @@ namespace Opm { } const double dt = simulator_.timeStepSize(); const Group& fieldGroup = this->schedule().getGroup("FIELD", reportStepIdx); - WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup, - this->schedule_, - regionalAveragePressureCalculator_, - reportStepIdx, - dt, - this->wellState(), - this->groupState()); + this->wgHelper().updateGpMaintTargetForGroups(fieldGroup, + regionalAveragePressureCalculator_, + dt, + this->groupState()); } this->updateAndCommunicateGroupData(reportStepIdx, @@ -513,9 +505,11 @@ namespace Opm { try { well->scaleSegmentRatesAndPressure(this->wellState()); well->calculateExplicitQuantities(simulator_, this->wellState(), local_deferredLogger); - well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), local_deferredLogger); + well->updateWellStateWithTarget(simulator_, this->wgHelper(), this->wellState(), local_deferredLogger); well->updatePrimaryVariables(simulator_, this->wellState(), local_deferredLogger); - well->solveWellEquation(simulator_, this->wellState(), this->groupState(), local_deferredLogger); + well->solveWellEquation( + simulator_, this->wgHelper(), this->wellState(), local_deferredLogger + ); } catch (const std::exception& e) { const std::string msg = "Compute initial well solution for new well " + well->name() + " failed. Continue with zero initial rates"; local_deferredLogger.warning("WELL_INITIAL_SOLVE_FAILED", msg); @@ -554,11 +548,10 @@ namespace Opm { Scalar well_efficiency_factor = wellEcl.getEfficiencyFactor() * this->wellState().getGlobalEfficiencyScalingFactor(well_name); - WellGroupHelpersType::accumulateGroupEfficiencyFactor(this->schedule().getGroup(wellEcl.groupName(), - timeStepIdx), - this->schedule(), - timeStepIdx, - well_efficiency_factor); + this->wgHelper().accumulateGroupEfficiencyFactor( + this->schedule().getGroup(wellEcl.groupName(), timeStepIdx), + well_efficiency_factor + ); well->setWellEfficiencyFactor(well_efficiency_factor); well->setVFPProperties(this->vfp_properties_.get()); @@ -590,8 +583,8 @@ namespace Opm { gaslift_.initGliftEclWellMap(well_container_, ecl_well_map); well->wellTesting(simulator_, simulationTime, + this->wgHelper(), this->wellState(), - this->groupState(), this->wellTestState(), ecl_well_map, this->well_open_times_, @@ -1366,8 +1359,7 @@ namespace Opm { auto target = WellGroupControls::getAutoChokeGroupProductionTargetRate( group.name(), parentGroup, - well_state, - group_state, + this->wgHelper(), this->schedule(), summary_state, resv_coeff, @@ -1402,7 +1394,9 @@ namespace Opm { const Well& well_ecl = this->wells_ecl_[well->indexOfWell()]; const auto inj_controls = Well::InjectionControls(0); const auto prod_controls = well_ecl.productionControls(summary_state); - well->iterateWellEqWithSwitching(this->simulator_, dt, inj_controls, prod_controls, well_state, group_state, local_deferredLogger, false, false); + well->iterateWellEqWithSwitching( + this->simulator_, dt, inj_controls, prod_controls, this->wgHelper(), this->wellState(), local_deferredLogger, false, false + ); rate = -tcalc.calcModeRateFromRates(ws.surface_rates); group_rate += rate; } @@ -1490,7 +1484,9 @@ namespace Opm { 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->wellState(), this->groupState(), local_deferredLogger); + well->prepareWellBeforeAssembling( + this->simulator_, dt, this->wgHelper(), this->wellState(), local_deferredLogger + ); } } @@ -1512,7 +1508,7 @@ namespace Opm { { OPM_TIMEFUNCTION(); for (auto& well : well_container_) { - well->assembleWellEq(simulator_, dt, this->wellState(), this->groupState(), deferred_logger); + well->assembleWellEq(simulator_, dt, this->wgHelper(), deferred_logger); } } @@ -1524,7 +1520,9 @@ namespace Opm { { OPM_TIMEFUNCTION(); for (auto& well : well_container_) { - well->prepareWellBeforeAssembling(simulator_, dt, this->wellState(), this->groupState(), deferred_logger); + well->prepareWellBeforeAssembling( + simulator_, dt, this->wgHelper(), this->wellState(), deferred_logger + ); } } @@ -1806,7 +1804,9 @@ namespace Opm { OPM_BEGIN_PARALLEL_TRY_CATCH() for (const auto& well : well_container_) { const auto mode = WellInterface::IndividualOrGroup::Group; - const bool changed_well = well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger); + const bool changed_well = well->updateWellControl( + simulator_, mode, this->wgHelper(), this->wellState(), deferred_logger + ); if (changed_well) { changed_well_to_group = changed_well || changed_well_to_group; } @@ -1829,7 +1829,9 @@ namespace Opm { OPM_BEGIN_PARALLEL_TRY_CATCH() for (const auto& well : well_container_) { const auto mode = WellInterface::IndividualOrGroup::Individual; - const bool changed_well = well->updateWellControl(simulator_, mode, this->wellState(), this->groupState(), deferred_logger); + const bool changed_well = well->updateWellControl( + simulator_, mode, this->wgHelper(), this->wellState(), deferred_logger + ); if (changed_well) { changed_well_individual = changed_well || changed_well_individual; } @@ -1901,7 +1903,9 @@ namespace Opm { 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->wellState(), this->groupState(), deferred_logger); + well->prepareWellBeforeAssembling( + this->simulator_, dt, this->wgHelper(), this->wellState(), deferred_logger + ); } } } @@ -1937,8 +1941,9 @@ namespace Opm { if (ws.production_cmode == Well::ProducerCMode::GRUP || ws.injection_cmode == Well::InjectorCMode::GRUP) { - well->updateWellStateWithTarget(simulator_, this->groupState(), - this->wellState(), deferred_logger); + well->updateWellStateWithTarget( + simulator_, this->wgHelper(), this->wellState(), deferred_logger + ); } } OPM_END_PARALLEL_TRY_CATCH("BlackoilWellModel::updateAndCommunicate failed: ", @@ -2007,6 +2012,7 @@ namespace Opm { const auto wasClosed = wellTestState.well_is_closed(wname); well->checkWellOperability(simulator_, this->wellState(), + this->wgHelper(), local_deferredLogger); const bool under_zero_target = well->wellUnderZeroGroupRateTarget(this->simulator_, @@ -2091,7 +2097,7 @@ namespace Opm { std::string cur_exc_msg; auto cur_exc_type = ExceptionType::NONE; try { - well->computeWellPotentials(simulator_, well_state_copy, potentials, deferred_logger); + well->computeWellPotentials(simulator_, well_state_copy, this->wgHelper(), potentials, deferred_logger); } // catch all possible exception and store type and message. OPM_PARALLEL_CATCH_CLAUSE(cur_exc_type, cur_exc_msg); @@ -2185,7 +2191,9 @@ namespace Opm { for (const auto& well : well_container_) { auto& events = this->wellState().well(well->indexOfWell()).events; if (events.hasEvent(WellState::event_mask)) { - well->updateWellStateWithTarget(simulator_, this->groupState(), this->wellState(), deferred_logger); + well->updateWellStateWithTarget( + simulator_, this->wgHelper(), this->wellState(), deferred_logger + ); well->updatePrimaryVariables(simulator_, this->wellState(), deferred_logger); // There is no new well control change input within a report step, // so next time step, the well does not consider to have effective events anymore. @@ -2198,7 +2206,9 @@ namespace Opm { // solve the well equation initially to improve the initial solution of the well model if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) { try { - well->solveWellEquation(simulator_, this->wellState(), this->groupState(), deferred_logger); + well->solveWellEquation( + simulator_, this->wgHelper(), this->wellState(), deferred_logger + ); } catch (const std::exception& e) { const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the previous rates"; deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg); diff --git a/opm/simulators/wells/FractionCalculator.cpp b/opm/simulators/wells/FractionCalculator.cpp index 0dcf358c64..f589ae97fe 100644 --- a/opm/simulators/wells/FractionCalculator.cpp +++ b/opm/simulators/wells/FractionCalculator.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include #include #include @@ -37,8 +37,7 @@ namespace Opm::WGHelpers { template FractionCalculator:: FractionCalculator(const Schedule& schedule, - const WellState& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, const int report_step, const GuideRate* guide_rate, @@ -46,8 +45,7 @@ FractionCalculator(const Schedule& schedule, const bool is_producer, const Phase injection_phase) : schedule_(schedule) - , well_state_(well_state) - , group_state_(group_state) + , wgHelper_(wgHelper) , summary_state_(summary_state) , report_step_(report_step) , guide_rate_(guide_rate) @@ -129,11 +127,11 @@ guideRateSum(const Group& group, for (const std::string& child_group : group.groups()) { bool included = (child_group == always_included_child); if (is_producer_) { - const auto ctrl = this->group_state_.production_control(child_group); + const auto ctrl = this->wgHelper().groupState().production_control(child_group); included |= (ctrl == Group::ProductionCMode::FLD) || (ctrl == Group::ProductionCMode::NONE); } else { - const auto ctrl = this->group_state_.injection_control(child_group, + const auto ctrl = this->wgHelper().groupState().injection_control(child_group, this->injection_phase_); included |= (ctrl == Group::InjectionCMode::FLD) || (ctrl == Group::InjectionCMode::NONE); @@ -149,9 +147,9 @@ guideRateSum(const Group& group, for (const std::string& child_well : group.wells()) { bool included = (child_well == always_included_child); if (is_producer_) { - included |= well_state_.isProductionGrup(child_well); + included |= this->wgHelper().wellState().isProductionGrup(child_well); } else { - included |= well_state_.isInjectionGrup(child_well); + included |= this->wgHelper().wellState().isInjectionGrup(child_well); } if (included) { number_of_included_well_or_groups++; @@ -169,8 +167,7 @@ guideRate(const std::string& name, const bool always_use_potentials) { if (schedule_.hasWell(name, report_step_)) { - return WellGroupHelpers::getGuideRate(name, schedule_, well_state_, group_state_, - report_step_, guide_rate_, target_); + return this->wgHelper().getGuideRate(name, target_); } else { if (groupControlledWells(name, always_included_child) > 0) { if (is_producer_ && guide_rate_->has(name) && !always_use_potentials) { @@ -196,14 +193,9 @@ int FractionCalculator:: groupControlledWells(const std::string& group_name, const std::string& always_included_child) { - return WellGroupHelpers::groupControlledWells(schedule_, - well_state_, - this->group_state_, - report_step_, - group_name, - always_included_child, - is_producer_, - injection_phase_); + return this->wgHelper().groupControlledWells( + group_name, always_included_child, is_producer_, injection_phase_ + ); } template @@ -211,9 +203,7 @@ GuideRate::RateVector FractionCalculator:: getGroupRateVector(const std::string& group_name) { assert(is_producer_); - return WellGroupHelpers::getProductionGroupRateVector(this->group_state_, - this->well_state_.phaseUsageInfo(), - group_name); + return this->wgHelper().getProductionGroupRateVector(group_name); } template class FractionCalculator; diff --git a/opm/simulators/wells/FractionCalculator.hpp b/opm/simulators/wells/FractionCalculator.hpp index 941dba81c2..49a06474b4 100644 --- a/opm/simulators/wells/FractionCalculator.hpp +++ b/opm/simulators/wells/FractionCalculator.hpp @@ -28,6 +28,7 @@ namespace Opm { template class GroupState; class Schedule; +template class WellGroupHelper; template class WellState; } @@ -37,9 +38,10 @@ template class FractionCalculator { public: + using WellGroupHelperType = WellGroupHelper; + FractionCalculator(const Schedule& schedule, - const WellState& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, const int report_step, const GuideRate* guide_rate, @@ -66,9 +68,10 @@ private: int groupControlledWells(const std::string& group_name, const std::string& always_included_child); GuideRate::RateVector getGroupRateVector(const std::string& group_name); + const WellGroupHelperType& wgHelper() const { return wgHelper_; } + const Schedule& schedule_; - const WellState& well_state_; - const GroupState& group_state_; + const WellGroupHelperType& wgHelper_; const SummaryState& summary_state_; int report_step_; const GuideRate* guide_rate_; diff --git a/opm/simulators/wells/GasLiftSingleWell_impl.hpp b/opm/simulators/wells/GasLiftSingleWell_impl.hpp index 60aa414023..0c8f89f304 100644 --- a/opm/simulators/wells/GasLiftSingleWell_impl.hpp +++ b/opm/simulators/wells/GasLiftSingleWell_impl.hpp @@ -146,8 +146,10 @@ GasLiftSingleWell:: computeBhpAtThpLimit_(Scalar alq, bool debug_output) const { OPM_TIMEFUNCTION(); + const auto& wgHelper = this->simulator_.problem().wellModel().wgHelper(); auto bhp_at_thp_limit = this->well_.computeBhpAtThpLimitProdWithAlq( this->simulator_, + wgHelper, this->summary_state_, alq, this->deferred_logger_, diff --git a/opm/simulators/wells/GroupEconomicLimitsChecker.cpp b/opm/simulators/wells/GroupEconomicLimitsChecker.cpp index 64d5550b37..1c56bd008f 100644 --- a/opm/simulators/wells/GroupEconomicLimitsChecker.cpp +++ b/opm/simulators/wells/GroupEconomicLimitsChecker.cpp @@ -31,7 +31,7 @@ #include #include -#include +#include #include @@ -79,13 +79,9 @@ GroupEconomicLimitsChecker(const BlackoilWellModelGeneric& auto phase_idx = this->phase_idx_map_[i]; this->phase_idx_reverse_map_[phase_idx] = static_cast(i); auto phase_pos = this->well_model_.phaseUsage().canonicalToActivePhaseIdx(phase_idx); - Scalar production_rate = WellGroupHelpers::sumWellSurfaceRates(this->group_, - this->schedule_, - this->well_state_, - this->report_step_idx_, - phase_pos, - /*isInjector*/false, - well_model_.summaryState()); + Scalar production_rate = this->well_model_.wgHelper().sumWellSurfaceRates( + this->group_, phase_pos, /*isInjector*/false + ); this->production_rates_[i] = this->well_model_.comm().sum(production_rate); } } diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 122af65b8e..53a54717ab 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -52,6 +52,7 @@ namespace Opm { using typename Base::SparseMatrixAdapter; using typename Base::FluidState; using typename Base::WellStateType; + using typename Base::WellGroupHelperType; using Base::has_solvent; using Base::has_polymer; @@ -90,7 +91,7 @@ namespace Opm { /// updating the well state based the current control mode void updateWellStateWithTarget(const Simulator& simulator, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger) const override; @@ -119,6 +120,7 @@ namespace Opm { /// computing the well potentials for group control void computeWellPotentials(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) override; @@ -160,6 +162,7 @@ namespace Opm { std::optional computeBhpAtThpLimitProdWithAlq(const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, const Scalar alq_value, DeferredLogger& deferred_logger, @@ -234,6 +237,7 @@ namespace Opm { DeferredLogger& deferred_logger) const; void computeWellRatesAtBhpLimit(const Simulator& simulator, + const WellGroupHelperType& wgHelper, std::vector& well_flux, DeferredLogger& deferred_logger) const; @@ -244,16 +248,18 @@ namespace Opm { void computeWellRatesWithBhpIterations(const Simulator& simulator, const Scalar& bhp, + const WellGroupHelperType& wgHelper, std::vector& well_flux, DeferredLogger& deferred_logger) const override; std::vector computeWellPotentialWithTHP(const WellStateType& well_state, const Simulator& simulator, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) const; bool computeWellPotentialsImplicit(const Simulator& simulator, - const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) const; @@ -263,16 +269,16 @@ namespace Opm { const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) override; bool iterateWellEqWithSwitching(const Simulator& simulator, const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger, const bool fixed_control = false, const bool fixed_status = false) override; @@ -305,11 +311,13 @@ namespace Opm { std::optional computeBhpAtThpLimitProd(const WellStateType& well_state, const Simulator& ebos_simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const; std::optional computeBhpAtThpLimitInj(const Simulator& ebos_simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const; @@ -323,6 +331,7 @@ namespace Opm { // check whether the well is operable under THP limit with current reservoir condition void checkOperabilityUnderTHPLimit(const Simulator& ebos_simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) override; // updating the inflow based on the current reservoir condition diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index d7c391adad..fdcc286879 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -182,11 +182,11 @@ namespace Opm void MultisegmentWell:: updateWellStateWithTarget(const Simulator& simulator, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger) const { - Base::updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger); + Base::updateWellStateWithTarget(simulator, wgHelper, well_state, deferred_logger); // scale segment rates based on the wellRates // and segment pressure based on bhp this->scaleSegmentRatesWithWellRates(this->segments_.inlets(), @@ -296,6 +296,7 @@ namespace Opm MultisegmentWell:: computeWellPotentials(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) { @@ -309,7 +310,7 @@ namespace Opm debug_cost_counter_ = 0; bool converged_implicit = false; if (this->param_.local_well_solver_control_switching_) { - converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger); + converged_implicit = computeWellPotentialsImplicit(simulator, wgHelper, well_potentials, deferred_logger); if (!converged_implicit) { deferred_logger.debug("Implicit potential calculations failed for well " + this->name() + ", reverting to original aproach."); @@ -319,10 +320,10 @@ namespace Opm // does the well have a THP related constraint? const auto& summaryState = simulator.vanguard().summaryState(); if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) { - computeWellRatesAtBhpLimit(simulator, well_potentials, deferred_logger); + computeWellRatesAtBhpLimit(simulator, wgHelper, well_potentials, deferred_logger); } else { well_potentials = computeWellPotentialWithTHP( - well_state, simulator, deferred_logger); + well_state, simulator, wgHelper, deferred_logger); } } deferred_logger.debug("Cost in iterations of finding well potential for well " @@ -340,15 +341,16 @@ namespace Opm void MultisegmentWell:: computeWellRatesAtBhpLimit(const Simulator& simulator, + const WellGroupHelperType& wgHelper, std::vector& well_flux, DeferredLogger& deferred_logger) const { if (this->well_ecl_.isInjector()) { const auto controls = this->well_ecl_.injectionControls(simulator.vanguard().summaryState()); - computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger); + computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, wgHelper, well_flux, deferred_logger); } else { const auto controls = this->well_ecl_.productionControls(simulator.vanguard().summaryState()); - computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, well_flux, deferred_logger); + computeWellRatesWithBhpIterations(simulator, controls.bhp_limit, wgHelper, well_flux, deferred_logger); } } @@ -404,6 +406,7 @@ namespace Opm MultisegmentWell:: computeWellRatesWithBhpIterations(const Simulator& simulator, const Scalar& bhp, + const WellGroupHelperType& wgHelper, std::vector& well_flux, DeferredLogger& deferred_logger) const { @@ -415,9 +418,10 @@ namespace Opm well_copy.debug_cost_counter_ = 0; + WellGroupHelperType wgHelper_copy = wgHelper; // store a copy of the well state, we don't want to update the real well state - WellStateType well_state_copy = simulator.problem().wellModel().wellState(); - const auto& group_state = simulator.problem().wellModel().groupState(); + WellStateType well_state_copy = wgHelper_copy.wellState(); + auto guard = wgHelper_copy.pushWellState(well_state_copy); auto& ws = well_state_copy.well(this->index_of_well_); // Get the current controls. @@ -459,8 +463,8 @@ namespace Opm well_copy.calculateExplicitQuantities(simulator, well_state_copy, deferred_logger); const double dt = simulator.timeStepSize(); // iterate to get a solution at the given bhp. - well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, - deferred_logger); + well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, wgHelper_copy, + well_state_copy, deferred_logger); // compute the potential and store in the flux vector. well_flux.clear(); @@ -479,6 +483,7 @@ namespace Opm MultisegmentWell:: computeWellPotentialWithTHP(const WellStateType& well_state, const Simulator& simulator, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) const { std::vector potentials(this->number_of_phases_, 0.0); @@ -486,12 +491,12 @@ namespace Opm const auto& well = this->well_ecl_; if (well.isInjector()) { - auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger); + auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, wgHelper, summary_state, deferred_logger); if (bhp_at_thp_limit) { const auto& controls = well.injectionControls(summary_state); const Scalar bhp = std::min(*bhp_at_thp_limit, static_cast(controls.bhp_limit)); - computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, potentials, deferred_logger); deferred_logger.debug("Converged thp based potential calculation for well " + this->name() + ", at bhp = " + std::to_string(bhp)); } else { @@ -500,16 +505,16 @@ namespace Opm + this->name() + ". Instead the bhp based value is used"); const auto& controls = well.injectionControls(summary_state); const Scalar bhp = controls.bhp_limit; - computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, potentials, deferred_logger); } } else { auto bhp_at_thp_limit = computeBhpAtThpLimitProd( - well_state, simulator, summary_state, deferred_logger); + well_state, simulator, wgHelper, summary_state, deferred_logger); if (bhp_at_thp_limit) { const auto& controls = well.productionControls(summary_state); const Scalar bhp = std::max(*bhp_at_thp_limit, static_cast(controls.bhp_limit)); - computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, potentials, deferred_logger); deferred_logger.debug("Converged thp based potential calculation for well " + this->name() + ", at bhp = " + std::to_string(bhp)); } else { @@ -518,7 +523,7 @@ namespace Opm + this->name() + ". Instead the bhp based value is used"); const auto& controls = well.productionControls(summary_state); const Scalar bhp = controls.bhp_limit; - computeWellRatesWithBhpIterations(simulator, bhp, potentials, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, potentials, deferred_logger); } } @@ -529,7 +534,7 @@ namespace Opm bool MultisegmentWell:: computeWellPotentialsImplicit(const Simulator& simulator, - const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) const { @@ -539,9 +544,10 @@ namespace Opm MultisegmentWell well_copy(*this); well_copy.debug_cost_counter_ = 0; + WellGroupHelperType wgHelper_copy = wgHelper; // store a copy of the well state, we don't want to update the real well state - WellStateType well_state_copy = well_state; - const auto& group_state = simulator.problem().wellModel().groupState(); + WellStateType well_state_copy = wgHelper_copy.wellState(); + auto guard = wgHelper_copy.pushWellState(well_state_copy); auto& ws = well_state_copy.well(this->index_of_well_); // get current controls @@ -579,9 +585,13 @@ namespace Opm // solve equations bool converged = false; if (this->well_ecl_.isProducer()) { - converged = well_copy.solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger); + converged = well_copy.solveWellWithOperabilityCheck( + simulator, dt, inj_controls, prod_controls, wgHelper_copy, well_state_copy, deferred_logger + ); } else { - converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger); + converged = well_copy.iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper_copy, well_state_copy, deferred_logger + ); } // fetch potentials (sign is updated on the outside). @@ -1432,13 +1442,14 @@ namespace Opm MultisegmentWell:: checkOperabilityUnderTHPLimit(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) { const auto& summaryState = simulator.vanguard().summaryState(); const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd( - well_state, simulator, summaryState, deferred_logger) - : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger); + well_state, simulator, wgHelper, summaryState, deferred_logger) + : computeBhpAtThpLimitInj(simulator, wgHelper, summaryState, deferred_logger); if (obtain_bhp) { this->operability_status_.can_obtain_bhp_with_thp_limit = true; @@ -1486,10 +1497,11 @@ namespace Opm const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { + const auto& group_state = wgHelper.groupState(); if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true; const int max_iter_number = this->param_.max_inner_iter_ms_wells_; @@ -1610,12 +1622,13 @@ namespace Opm const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger, const bool fixed_control /*false*/, const bool fixed_status /*false*/) { + const auto& group_state = wgHelper.groupState(); const int max_iter_number = this->param_.max_inner_iter_ms_wells_; { @@ -1664,10 +1677,10 @@ namespace Opm ++its_since_last_switch; if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){ const Scalar wqTotal = this->primary_variables_.getWQTotal().value(); - bool changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state, - inj_controls, prod_controls, wqTotal, - deferred_logger, fixed_control, - fixed_status); + bool changed = this->updateWellControlAndStatusLocalIteration( + simulator, wgHelper, inj_controls, prod_controls, wqTotal, + well_state, deferred_logger, fixed_control, fixed_status + ); if (changed) { its_since_last_switch = 0; ++switch_count; @@ -2059,11 +2072,13 @@ namespace Opm MultisegmentWell:: computeBhpAtThpLimitProd(const WellStateType& well_state, const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const { return this->MultisegmentWell::computeBhpAtThpLimitProdWithAlq( simulator, + wgHelper, summary_state, this->getALQ(well_state), deferred_logger, @@ -2076,6 +2091,7 @@ namespace Opm std::optional::Scalar> MultisegmentWell:: computeBhpAtThpLimitProdWithAlq(const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, const Scalar alq_value, DeferredLogger& deferred_logger, @@ -2109,12 +2125,12 @@ namespace Opm if (!iterate_if_no_solution) return std::nullopt; - auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) { + auto fratesIter = [this, &simulator, &wgHelper, &deferred_logger](const Scalar bhp) { // Solver the well iterations to see if we are // able to get a solution with an update // solution std::vector rates(3); - computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, rates, deferred_logger); return rates; }; @@ -2132,6 +2148,7 @@ namespace Opm std::optional::Scalar> MultisegmentWell:: computeBhpAtThpLimitInj(const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const { @@ -2159,12 +2176,12 @@ namespace Opm if (bhpAtLimit) return bhpAtLimit; - auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) { + auto fratesIter = [this, &simulator, &wgHelper, &deferred_logger](const Scalar bhp) { // Solver the well iterations to see if we are // able to get a solution with an update // solution std::vector rates(3); - computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, rates, deferred_logger); return rates; }; diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index c01efe54ce..1bb974661b 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -77,6 +77,7 @@ namespace Opm using typename Base::SparseMatrixAdapter; using typename Base::FluidState; using typename Base::RateVector; + using typename Base::WellGroupHelperType; using Base::has_solvent; using Base::has_zFraction; @@ -161,6 +162,7 @@ namespace Opm /// computing the well potentials for group control void computeWellPotentials(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) /* const */ override; @@ -197,8 +199,8 @@ namespace Opm const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) override; // iterate well equations including control switching @@ -206,20 +208,22 @@ namespace Opm const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger, const bool fixed_control = false, const bool fixed_status = false) override; /* returns BHP */ Scalar computeWellRatesAndBhpWithThpAlqProd(const Simulator& ebos_simulator, + const WellGroupHelperType& wgHelper, const SummaryState &summary_state, DeferredLogger& deferred_logger, std::vector& potentials, Scalar alq) const; void computeWellRatesWithThpAlqProd(const Simulator& ebos_simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger, std::vector& potentials, @@ -227,6 +231,7 @@ namespace Opm std::optional computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, const Scalar alq_value, DeferredLogger& deferred_logger, @@ -309,16 +314,18 @@ namespace Opm void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator, const Scalar& bhp, + const WellGroupHelperType& wgHelper, std::vector& well_flux, DeferredLogger& deferred_logger) const override; std::vector computeWellPotentialWithTHP(const Simulator& ebosSimulator, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger, const WellStateType& well_state) const; bool computeWellPotentialsImplicit(const Simulator& ebos_simulator, - const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) const; @@ -379,6 +386,7 @@ namespace Opm // check whether the well is operable under THP limit with current reservoir condition void checkOperabilityUnderTHPLimit(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) override; // updating the inflow based on the current reservoir condition @@ -445,11 +453,13 @@ namespace Opm std::optional computeBhpAtThpLimitProd(const WellStateType& well_state, const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const; std::optional computeBhpAtThpLimitInj(const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 724ad60722..ff0b5d4a68 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1036,11 +1036,12 @@ namespace Opm StandardWell:: checkOperabilityUnderTHPLimit(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) { const auto& summaryState = simulator.vanguard().summaryState(); - const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, summaryState, deferred_logger) - : computeBhpAtThpLimitInj(simulator, summaryState, deferred_logger); + const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(well_state, simulator, wgHelper, summaryState, deferred_logger) + : computeBhpAtThpLimitInj(simulator, wgHelper, summaryState, deferred_logger); if (obtain_bhp) { this->operability_status_.can_obtain_bhp_with_thp_limit = true; @@ -1516,6 +1517,7 @@ namespace Opm StandardWell:: computeWellRatesWithBhpIterations(const Simulator& simulator, const Scalar& bhp, + const WellGroupHelperType& wgHelper, std::vector& well_flux, DeferredLogger& deferred_logger) const { @@ -1527,8 +1529,12 @@ namespace Opm // iterate to get a more accurate well density // create a copy of the well_state to use. If the operability checking is sucessful, we use this one // to replace the original one - WellStateType well_state_copy = simulator.problem().wellModel().wellState(); - const auto& group_state = simulator.problem().wellModel().groupState(); + WellGroupHelperType wgHelper_copy = wgHelper; + WellStateType well_state_copy = wgHelper_copy.wellState(); + // Ensure that wgHelper_copy uses well_state_copy as WellState for the rest of this function, + // and the guard ensures that the original well state is restored at scope exit, i.e. at + // the end of this function. + auto guard = wgHelper_copy.pushWellState(well_state_copy); // Get the current controls. const auto& summary_state = simulator.vanguard().summaryState(); @@ -1561,7 +1567,9 @@ namespace Opm well_copy.computeAccumWell(); const double dt = simulator.timeStepSize(); - const bool converged = well_copy.iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger); + const bool converged = well_copy.iterateWellEqWithControl( + simulator, dt, inj_controls, prod_controls, wgHelper_copy, well_state_copy, deferred_logger + ); if (!converged) { const std::string msg = " well " + name() + " did not get converged during well potential calculations " " potentials are computed based on unconverged solution"; @@ -1579,6 +1587,7 @@ namespace Opm std::vector::Scalar> StandardWell:: computeWellPotentialWithTHP(const Simulator& simulator, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger, const WellStateType& well_state) const { @@ -1588,7 +1597,7 @@ namespace Opm const auto& well = this->well_ecl_; if (well.isInjector()){ const auto& controls = this->well_ecl_.injectionControls(summary_state); - auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, summary_state, deferred_logger); + auto bhp_at_thp_limit = computeBhpAtThpLimitInj(simulator, wgHelper, summary_state, deferred_logger); if (bhp_at_thp_limit) { const Scalar bhp = std::min(*bhp_at_thp_limit, static_cast(controls.bhp_limit)); @@ -1602,7 +1611,7 @@ namespace Opm } } else { computeWellRatesWithThpAlqProd( - simulator, summary_state, + simulator, wgHelper, summary_state, deferred_logger, potentials, this->getALQ(well_state) ); } @@ -1614,7 +1623,7 @@ namespace Opm bool StandardWell:: computeWellPotentialsImplicit(const Simulator& simulator, - const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) const { @@ -1624,8 +1633,12 @@ namespace Opm StandardWell well_copy(*this); // store a copy of the well state, we don't want to update the real well state - WellStateType well_state_copy = well_state; - const auto& group_state = simulator.problem().wellModel().groupState(); + WellStateType well_state_copy = wgHelper.wellState(); + WellGroupHelperType wgHelper_copy = wgHelper; + // Ensure that wgHelper_copy uses well_state_copy as WellState for the rest of this function, + // and the guard ensures that the original well state is restored at scope exit, i.e. at + // the end of this function. + auto guard = wgHelper_copy.pushWellState(well_state_copy); auto& ws = well_state_copy.well(this->index_of_well_); // get current controls @@ -1663,9 +1676,13 @@ namespace Opm // iterate to get a solution at the given bhp. bool converged = false; if (this->well_ecl_.isProducer()) { - converged = well_copy.solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger); + converged = well_copy.solveWellWithOperabilityCheck( + simulator, dt, inj_controls, prod_controls, wgHelper_copy, well_state_copy, deferred_logger + ); } else { - converged = well_copy.iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state_copy, group_state, deferred_logger); + converged = well_copy.iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper_copy, well_state_copy, deferred_logger + ); } // fetch potentials (sign is updated on the outside). @@ -1693,6 +1710,7 @@ namespace Opm typename StandardWell::Scalar StandardWell:: computeWellRatesAndBhpWithThpAlqProd(const Simulator &simulator, + const WellGroupHelperType& wgHelper, const SummaryState &summary_state, DeferredLogger& deferred_logger, std::vector& potentials, @@ -1700,7 +1718,7 @@ namespace Opm { Scalar bhp; auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq( - simulator, summary_state, alq, deferred_logger, /*iterate_if_no_solution */ true); + simulator, wgHelper, summary_state, alq, deferred_logger, /*iterate_if_no_solution */ true); if (bhp_at_thp_limit) { const auto& controls = this->well_ecl_.productionControls(summary_state); bhp = std::max(*bhp_at_thp_limit, @@ -1722,6 +1740,7 @@ namespace Opm void StandardWell:: computeWellRatesWithThpAlqProd(const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger, std::vector& potentials, @@ -1729,6 +1748,7 @@ namespace Opm { /*double bhp =*/ computeWellRatesAndBhpWithThpAlqProd(simulator, + wgHelper, summary_state, deferred_logger, potentials, @@ -1740,6 +1760,7 @@ namespace Opm StandardWell:: computeWellPotentials(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) // const { @@ -1755,7 +1776,9 @@ namespace Opm // group controlled wells with defaulted guiderates will have zero targets as // the potentials are used to compute the well fractions. if (this->param_.local_well_solver_control_switching_ && !(this->changed_to_open_this_step_ && this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger))) { - converged_implicit = computeWellPotentialsImplicit(simulator, well_state, well_potentials, deferred_logger); + converged_implicit = computeWellPotentialsImplicit( + simulator, wgHelper, well_potentials, deferred_logger + ); } if (!converged_implicit) { // does the well have a THP related constraint? @@ -1777,10 +1800,10 @@ namespace Opm bhp = std::min(ws.bhp, bhp); assert(std::abs(bhp) != std::numeric_limits::max()); - computeWellRatesWithBhpIterations(simulator, bhp, well_potentials, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, well_potentials, deferred_logger); } else { // the well has a THP related constraint - well_potentials = computeWellPotentialWithTHP(simulator, deferred_logger, well_state); + well_potentials = computeWellPotentialWithTHP(simulator, wgHelper, deferred_logger, well_state); } } @@ -2210,10 +2233,12 @@ namespace Opm StandardWell:: computeBhpAtThpLimitProd(const WellStateType& well_state, const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const { return computeBhpAtThpLimitProdWithAlq(simulator, + wgHelper, summary_state, this->getALQ(well_state), deferred_logger, @@ -2224,6 +2249,7 @@ namespace Opm std::optional::Scalar> StandardWell:: computeBhpAtThpLimitProdWithAlq(const Simulator& simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, const Scalar alq_value, DeferredLogger& deferred_logger, @@ -2273,12 +2299,12 @@ namespace Opm if (!iterate_if_no_solution) return std::nullopt; - auto fratesIter = [this, &simulator, &deferred_logger](const Scalar bhp) { + auto fratesIter = [this, &simulator, &wgHelper, &deferred_logger](const Scalar bhp) { // Solver the well iterations to see if we are // able to get a solution with an update // solution std::vector rates(3); - computeWellRatesWithBhpIterations(simulator, bhp, rates, deferred_logger); + computeWellRatesWithBhpIterations(simulator, bhp, wgHelper, rates, deferred_logger); this->adaptRatesForVFP(rates); return rates; }; @@ -2310,9 +2336,15 @@ namespace Opm std::optional::Scalar> StandardWell:: computeBhpAtThpLimitInj(const Simulator& simulator, + [[maybe_unused]] const WellGroupHelperType& wgHelper, const SummaryState& summary_state, DeferredLogger& deferred_logger) const { + // Note: wgHelper parameter is currently unused in StandardWell but kept for consistency + // with MultisegmentWell::computeBhpAtThpLimitInj which uses it in fratesIter lambda. + // This maintains parallel API structure between well types and allows for future + // enhancements without breaking the interface. + // Make the frates() function. auto frates = [this, &simulator, &deferred_logger](const Scalar bhp) { // Not solving the well equations here, which means we are @@ -2345,10 +2377,11 @@ namespace Opm const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { + const auto& group_state = wgHelper.groupState(); updatePrimaryVariables(simulator, well_state, deferred_logger); const int max_iter = this->param_.max_inner_iter_wells_; @@ -2392,12 +2425,13 @@ namespace Opm const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger, const bool fixed_control /*false*/, const bool fixed_status /*false*/) { + const auto& group_state = wgHelper.groupState(); updatePrimaryVariables(simulator, well_state, deferred_logger); const int max_iter = this->param_.max_inner_iter_wells_; @@ -2436,9 +2470,10 @@ namespace Opm its_since_last_switch++; if (allow_switching && its_since_last_switch >= min_its_after_switch && status_switch_count < max_status_switch){ const Scalar wqTotal = this->primary_variables_.eval(WQTotal).value(); - changed = this->updateWellControlAndStatusLocalIteration(simulator, well_state, group_state, - inj_controls, prod_controls, wqTotal, - deferred_logger, fixed_control, fixed_status); + changed = this->updateWellControlAndStatusLocalIteration( + simulator, wgHelper, inj_controls, prod_controls, wqTotal, + well_state, deferred_logger, fixed_control, fixed_status + ); if (changed){ its_since_last_switch = 0; switch_count++; diff --git a/opm/simulators/wells/WellGroupConstraints.cpp b/opm/simulators/wells/WellGroupConstraints.cpp index e823f8dd3c..c0cc0dd160 100644 --- a/opm/simulators/wells/WellGroupConstraints.cpp +++ b/opm/simulators/wells/WellGroupConstraints.cpp @@ -27,7 +27,7 @@ #include -#include +#include #include #include @@ -38,15 +38,15 @@ template std::pair WellGroupConstraints:: checkGroupConstraintsInj(const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const Scalar efficiencyFactor, - const Schedule& schedule, const SummaryState& summaryState, const RateConvFunc& rateConverter, const bool check_guide_rate, DeferredLogger& deferred_logger) const { + const auto& well_state = wgHelper.wellState(); + // Translate injector type from control to Phase. const auto& well_controls = well_.wellEcl().injectionControls(summaryState); auto injectorType = well_controls.injector_type; @@ -77,65 +77,55 @@ checkGroupConstraintsInj(const Group& group, const auto& ws = well_state.well(well_.indexOfWell()); // Call check for the well's injection phase. - return WellGroupHelpers::checkGroupConstraintsInj(well_.name(), - well_.wellEcl().groupName(), - group, - well_state, - group_state, - well_.currentStep(), - well_.guideRate(), - ws.surface_rates.data(), - injectionPhase, - efficiencyFactor, - schedule, - summaryState, - resv_coeff, - check_guide_rate, - deferred_logger); + return wgHelper.checkGroupConstraintsInj( + well_.name(), + well_.wellEcl().groupName(), + group, + ws.surface_rates.data(), + injectionPhase, + efficiencyFactor, + resv_coeff, + check_guide_rate, + deferred_logger + ); } template std::pair WellGroupConstraints:: checkGroupConstraintsProd(const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, const RateConvFunc& rateConverter, const bool check_guide_rate, DeferredLogger& deferred_logger) const { + const auto& well_state = wgHelper.wellState(); + // Make conversion factors for RESV <-> surface rates. std::vector resv_coeff(well_.phaseUsage().numActivePhases(), 1.0); rateConverter(0, well_.pvtRegionIdx(), group.name(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS. const auto& ws = well_state.well(well_.indexOfWell()); - return WellGroupHelpers::checkGroupConstraintsProd(well_.name(), - well_.wellEcl().groupName(), - group, - well_state, - group_state, - well_.currentStep(), - well_.guideRate(), - ws.surface_rates.data(), - efficiencyFactor, - schedule, - summaryState, - resv_coeff, - check_guide_rate, - deferred_logger); + return wgHelper.checkGroupConstraintsProd(well_.name(), + well_.wellEcl().groupName(), + group, + ws.surface_rates.data(), + efficiencyFactor, + resv_coeff, + check_guide_rate, + deferred_logger + ); } template bool WellGroupConstraints:: -checkGroupConstraints(WellState& well_state, - const GroupState& group_state, +checkGroupConstraints(const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, const RateConvFunc& rateConverter, const bool check_guide_rate, + WellStateType& well_state, DeferredLogger& deferred_logger) const { const auto& well = well_.wellEcl(); @@ -156,9 +146,8 @@ checkGroupConstraints(WellState& well_state, const Scalar efficiencyFactor = well.getEfficiencyFactor() * well_state[well.name()].efficiency_scaling_factor; const std::pair group_constraint = - this->checkGroupConstraintsInj(group, well_state, - group_state, efficiencyFactor, - schedule, summaryState, + this->checkGroupConstraintsInj(group, wgHelper, efficiencyFactor, + summaryState, rateConverter, check_guide_rate, deferred_logger); @@ -189,10 +178,11 @@ checkGroupConstraints(WellState& well_state, const Scalar efficiencyFactor = well.getEfficiencyFactor() * well_state[well.name()].efficiency_scaling_factor; const std::pair group_constraint = - this->checkGroupConstraintsProd(group, well_state, - group_state, efficiencyFactor, - schedule, summaryState, - rateConverter, check_guide_rate, deferred_logger); + this->checkGroupConstraintsProd(group, wgHelper, + efficiencyFactor, + rateConverter, + check_guide_rate, + deferred_logger); // If a group constraint was broken, we set the current well control to // be GRUP. if (group_constraint.first) { diff --git a/opm/simulators/wells/WellGroupConstraints.hpp b/opm/simulators/wells/WellGroupConstraints.hpp index 5eaf171d78..85b8a3811c 100644 --- a/opm/simulators/wells/WellGroupConstraints.hpp +++ b/opm/simulators/wells/WellGroupConstraints.hpp @@ -40,6 +40,7 @@ enum class InjectorType; using RegionId = int; class Schedule; class SummaryState; +template class WellGroupHelper; template class WellInterfaceGeneric; template class WellState; @@ -54,22 +55,22 @@ public: const int, const std::optional&, std::vector&)>; + using WellGroupHelperType = WellGroupHelper; + using WellStateType = WellState; - bool checkGroupConstraints(WellState& well_state, - const GroupState& group_state, + bool checkGroupConstraints(const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, const RateConvFunc& rateConverter, const bool check_guide_rate, + WellStateType& well_state, DeferredLogger& deferred_logger) const; private: std::pair checkGroupConstraintsInj(const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const Scalar efficiencyFactor, - const Schedule& schedule, const SummaryState& summaryState, const RateConvFunc& rateConverter, const bool check_guide_rate, @@ -77,11 +78,8 @@ private: std::pair checkGroupConstraintsProd(const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, const RateConvFunc& rateConverter, const bool check_guide_rate, DeferredLogger& deferred_logger) const; diff --git a/opm/simulators/wells/WellGroupControls.cpp b/opm/simulators/wells/WellGroupControls.cpp index 0d45a97806..c084ecaf01 100644 --- a/opm/simulators/wells/WellGroupControls.cpp +++ b/opm/simulators/wells/WellGroupControls.cpp @@ -34,7 +34,7 @@ #include #include #include -#include +#include #include #include @@ -358,8 +358,7 @@ template std::pair WellGroupControls:: getAutoChokeGroupProductionTargetRate(const std::string& name, const Group& group, - const WellState& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, const std::vector& resv_coeff, @@ -368,6 +367,8 @@ getAutoChokeGroupProductionTargetRate(const std::string& name, const GuideRate* guideRate, DeferredLogger& deferred_logger) { + const auto& group_state = wgHelper.groupState(); + const auto& well_state = wgHelper.wellState(); const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name()); if (currentGroupControl == Group::ProductionCMode::FLD || currentGroupControl == Group::ProductionCMode::NONE) { @@ -377,7 +378,7 @@ getAutoChokeGroupProductionTargetRate(const std::string& name, // Produce share of parents control const auto& parent = schedule.getGroup(group.parent(), reportStepIdx); efficiencyFactor *= group.getGroupEfficiencyFactor(); - return getAutoChokeGroupProductionTargetRate(name, parent, well_state, group_state, + return getAutoChokeGroupProductionTargetRate(name, parent, wgHelper, schedule, summaryState, resv_coeff, efficiencyFactor, reportStepIdx, guideRate, deferred_logger); @@ -398,7 +399,7 @@ getAutoChokeGroupProductionTargetRate(const std::string& name, // gconsale may adjust the grat target. // the adjusted rates is send to the targetCalculator Scalar gratTargetFromSales = 0.0; - if (group_state.has_grat_sales_target(group.name())) + if (wgHelper.groupState().has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); WGHelpers::TargetCalculator tcalc(currentGroupControl, @@ -410,8 +411,7 @@ getAutoChokeGroupProductionTargetRate(const std::string& name, group.has_gpmaint_control(currentGroupControl)); WGHelpers::FractionCalculator fcalc(schedule, - well_state, - group_state, + wgHelper, summaryState, reportStepIdx, guideRate, @@ -433,8 +433,7 @@ getAutoChokeGroupProductionTargetRate(const std::string& name, ctrl = group.productionControls(summaryState); const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); - const auto chain = WellGroupHelpers::groupChainTopBot(name, group.name(), - schedule, reportStepIdx); + const auto chain = wgHelper.groupChainTopBot(name, group.name()); // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. const std::size_t num_ancestors = chain.size() - 1; Scalar target = orig_target; diff --git a/opm/simulators/wells/WellGroupControls.hpp b/opm/simulators/wells/WellGroupControls.hpp index 492830a8db..efb24541e8 100644 --- a/opm/simulators/wells/WellGroupControls.hpp +++ b/opm/simulators/wells/WellGroupControls.hpp @@ -40,6 +40,7 @@ enum class InjectorType; using RegionId = int; class Schedule; class SummaryState; +template class WellGroupHelper; template class WellInterfaceGeneric; template class WellState; @@ -52,6 +53,7 @@ public: using RateConvFunc = std::function&, std::vector&)>; + using WellGroupHelperType = WellGroupHelper; template void getGroupInjectionControl(const Group& group, @@ -101,9 +103,8 @@ public: DeferredLogger& deferred_logger) const; static std::pair getAutoChokeGroupProductionTargetRate(const std::string& name, - const Group& parent, - const WellState& well_state, - const GroupState& group_state, + const Group& group, + const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, const std::vector& resv_coeff, diff --git a/opm/simulators/wells/WellGroupHelper.cpp b/opm/simulators/wells/WellGroupHelper.cpp new file mode 100644 index 0000000000..e79ffe0526 --- /dev/null +++ b/opm/simulators/wells/WellGroupHelper.cpp @@ -0,0 +1,1964 @@ +/* + Copyright 2025 Equinor ASA + + 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 +WellGroupHelper::WellGroupHelper(WellState& well_state, + GroupState& group_state, + const Schedule& schedule, + const SummaryState& summary_state, + const GuideRate& guide_rate, + const PhaseUsageInfo& phase_usage_info) + : well_state_ {&well_state} + , group_state_ {&group_state} + , schedule_ {schedule} + , summary_state_ {summary_state} + , guide_rate_ {guide_rate} + , phase_usage_info_ {phase_usage_info} +{ +} + + +// --------------------------------------------------------------------- +// Public methods +// --------------------------------------------------------------------- + +template +void +WellGroupHelper::accumulateGroupEfficiencyFactor(const Group& group, + Scalar& factor) const +{ + factor *= group.getGroupEfficiencyFactor(); + if (group.parent() != "FIELD" && !group.parent().empty()) + this->accumulateGroupEfficiencyFactor(this->schedule_.getGroup(group.parent(), this->report_step_), + factor); +} + +template +std::pair +WellGroupHelper::checkGroupConstraintsInj(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + const Phase injection_phase, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + const bool check_guide_rate, + DeferredLogger& deferred_logger) const +{ + // When called for a well ('name' is a well name), 'parent' + // will be the name of 'group'. But if we recurse, 'name' and + // 'parent' will stay fixed while 'group' will be higher up + // in the group tree. + // efficiencyfactor is the well efficiency factor for the first group the well is + // part of. Later it is the accumulated factor including the group efficiency factor + // of the child of group. + + // NOTE: if name is a group, the rates-argument is the portion of the 'full' group rates that + // is potentially available for control by an ancestor, i.e. full rates minus reduction rates + + auto current_group_control = this->groupState().injection_control(group.name(), injection_phase); + + if (current_group_control == Group::InjectionCMode::FLD + || current_group_control == Group::InjectionCMode::NONE) { + // Return if we are not available for parent group. + if (!group.injectionGroupControlAvailable(injection_phase)) { + return std::make_pair(false, Scalar(1)); + } + // Otherwise: check production share of parent's control. + const auto& parent_group = this->schedule_.getGroup(group.parent(), this->report_step_); + return this->checkGroupConstraintsInj(name, + parent, + parent_group, + rates, + injection_phase, + efficiency_factor * group.getGroupEfficiencyFactor(), + resv_coeff, + check_guide_rate, + deferred_logger); + } + + // This can be false for FLD-controlled groups, we must therefore + // check for FLD first (done above). + if (!group.isInjectionGroup()) { + return std::make_pair(false, Scalar(1.0)); + } + + // If we are here, we are at the topmost group to be visited in the recursion. + // This is the group containing the control we will check against. + Scalar sales_target = 0; + if (this->schedule_[this->report_step_].gconsale().has(group.name())) { + const auto& gconsale + = this->schedule_[this->report_step_].gconsale().get(group.name(), this->summary_state_); + sales_target = gconsale.sales_target; + } + WGHelpers::InjectionTargetCalculator tcalc { + current_group_control, + this->phase_usage_info_, + resv_coeff, + group.name(), + sales_target, + this->groupState(), + injection_phase, + group.has_gpmaint_control(injection_phase, current_group_control), + deferred_logger}; + + WGHelpers::FractionCalculator fcalc {this->schedule_, + *this, + this->summary_state_, + this->report_step_, + &this->guide_rate_, + tcalc.guideTargetMode(), + /*is_production_group=*/false, + injection_phase}; + + auto local_fraction_lambda = [&](const std::string& child) { return fcalc.localFraction(child, name); }; + + auto local_reduction_lambda = [&](const std::string& group_name) { + const std::vector& group_target_reductions + = this->groupState().injection_reduction_rates(group_name); + return tcalc.calcModeRateFromRates(group_target_reductions); + }; + + auto local_current_rate_lambda = [&](const std::string& group_name) { + const std::vector& group_surface_rates + = this->groupState().injection_surface_rates(group_name); + return tcalc.calcModeRateFromRates(group_surface_rates); + }; + + const auto chain = this->groupChainTopBot(name, group.name()); + // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. + const std::size_t num_ancestors = chain.size() - 1; + // we need to find out the level where the current well is applied to the local reduction + std::size_t local_reduction_level = 0; + for (std::size_t ii = 1; ii < num_ancestors; ++ii) { + const int num_gr_ctrl = this->groupControlledWells(chain[ii], + /*always_included_child=*/"", + /*is_production_group=*/false, + injection_phase); + if (this->guide_rate_.has(chain[ii], injection_phase) && num_gr_ctrl > 0) { + local_reduction_level = ii; + } + } + + // check whether guide rate is violated + if (check_guide_rate) { + for (std::size_t ii = 1; ii < num_ancestors; ++ii) { + if (this->guide_rate_.has(chain[ii], injection_phase)) { + const auto& guided_group = chain[ii]; + const Scalar grefficiency + = this->schedule_.getGroup(guided_group, this->report_step_).getGroupEfficiencyFactor(); + const Scalar current_rate_fraction = grefficiency * local_current_rate_lambda(guided_group) + / local_current_rate_lambda(chain[ii - 1]); + const Scalar guiderate_fraction = local_fraction_lambda(guided_group); + // we add a factor here to avoid switching due to numerical instability + const Scalar factor = 1.01; + if (current_rate_fraction > (guiderate_fraction * factor)) { + return std::make_pair(true, guiderate_fraction / current_rate_fraction); + } + } + } + } + + if (this->schedule_.hasWell(name) + && this->wellState().well(name).group_target) { // for wells we already have computed the target + Scalar scale = 1.0; + const Scalar group_target_rate_available = *this->wellState().well(name).group_target; + const Scalar current_well_rate_available + = tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. + if (current_well_rate_available > 1e-12) { + scale = group_target_rate_available / current_well_rate_available; + } + return std::make_pair(current_well_rate_available > group_target_rate_available, scale); + } + + std::optional ctrl; + if (!group.has_gpmaint_control(injection_phase, current_group_control)) + ctrl = group.injectionControls(injection_phase, this->summary_state_); + + + const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); + // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. + // Then ... + // TODO finish explanation. + const Scalar current_rate_available + = tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. + + // Compute portion of target corresponding to current_rate_available + Scalar target = orig_target; + for (std::size_t ii = 0; ii < num_ancestors; ++ii) { + if ((ii == 0) || this->guide_rate_.has(chain[ii], injection_phase)) { + // Apply local reductions only at the control level + // (top) and for levels where we have a specified + // group guide rate. + if (local_reduction_level >= ii) { + target -= local_reduction_lambda(chain[ii]); + } + + // Add my reduction back at the level where it is included in the local reduction + if (local_reduction_level == ii) { + target += current_rate_available * efficiency_factor; + } + } + target *= local_fraction_lambda(chain[ii + 1]); + } + // Avoid negative target rates comming from too large local reductions. + const Scalar target_rate_available = std::max(Scalar(1e-12), target / efficiency_factor); + Scalar scale = 1.0; + if (current_rate_available > 1e-12) + scale = target_rate_available / current_rate_available; + + return std::make_pair(current_rate_available > target_rate_available, scale); +} + +template +std::pair +WellGroupHelper::checkGroupConstraintsProd(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + const bool check_guide_rate, + DeferredLogger& deferred_logger) const +{ + // When called for a well ('name' is a well name), 'parent' + // will be the name of 'group'. But if we recurse, 'name' and + // 'parent' will stay fixed while 'group' will be higher up + // in the group tree. + // efficiencyfactor is the well efficiency factor for the first group the well is + // part of. Later it is the accumulated factor including the group efficiency factor + // of the child of group. + + // NOTE: if name is a group, the rates-argument is the portion of the 'full' group rates that + // is potentially available for control by an ancestor, i.e. full rates minus reduction rates + OPM_TIMEFUNCTION(); + const Group::ProductionCMode& current_group_control = this->groupState().production_control(group.name()); + + if (current_group_control == Group::ProductionCMode::FLD + || current_group_control == Group::ProductionCMode::NONE) { + // Return if we are not available for parent group. + if (!group.productionGroupControlAvailable()) { + return std::make_pair(false, 1); + } + // Otherwise: check production share of parent's control. + const auto& parent_group = this->schedule_.getGroup(group.parent(), this->report_step_); + return this->checkGroupConstraintsProd(name, + parent, + parent_group, + rates, + efficiency_factor * group.getGroupEfficiencyFactor(), + resv_coeff, + check_guide_rate, + deferred_logger); + } + + // This can be false for FLD-controlled groups, we must therefore + // check for FLD first (done above). + if (!group.isProductionGroup()) { + return std::make_pair(false, Scalar(1.0)); + } + + // If we are here, we are at the topmost group to be visited in the recursion. + // This is the group containing the control we will check against. + + // gconsale may adjust the grat target. + // the adjusted rates is send to the targetCalculator + Scalar grat_target_from_sales = 0.0; + if (this->groupState().has_grat_sales_target(group.name())) + grat_target_from_sales = this->groupState().grat_sales_target(group.name()); + + WGHelpers::TargetCalculator tcalc {current_group_control, + this->phase_usage_info_, + resv_coeff, + grat_target_from_sales, + group.name(), + this->groupState(), + group.has_gpmaint_control(current_group_control)}; + + WGHelpers::FractionCalculator fcalc {this->schedule_, + *this, + this->summary_state_, + this->report_step_, + &this->guide_rate_, + tcalc.guideTargetMode(), + /*is_production_group=*/true, + /*injection_phase=*/Phase::OIL}; + + auto local_fraction_lambda = [&](const std::string& child) { return fcalc.localFraction(child, name); }; + + auto local_reduction_lambda = [&](const std::string& group_name) { + const std::vector& group_target_reductions + = this->groupState().production_reduction_rates(group_name); + return tcalc.calcModeRateFromRates(group_target_reductions); + }; + + auto local_current_rate_lambda = [&](const std::string& group_name) { + const std::vector& group_surface_rates = this->groupState().production_rates(group_name); + return tcalc.calcModeRateFromRates(group_surface_rates); + }; + + std::optional ctrl; + if (!group.has_gpmaint_control(current_group_control)) + ctrl = group.productionControls(this->summary_state_); + + const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); + // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. + // Then ... + // TODO finish explanation. + const Scalar current_rate_available + = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. + const auto chain = this->groupChainTopBot(name, group.name()); + // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. + const std::size_t num_ancestors = chain.size() - 1; + + // check whether guide rate is violated + if (check_guide_rate) { + for (std::size_t ii = 1; ii < num_ancestors; ++ii) { + if (this->guide_rate_.has(chain[ii])) { + const auto& guided_group = chain[ii]; + const Scalar grefficiency + = this->schedule_.getGroup(guided_group, this->report_step_).getGroupEfficiencyFactor(); + const Scalar current_rate_fraction = grefficiency * local_current_rate_lambda(guided_group) + / (local_current_rate_lambda(chain[ii - 1])); + const Scalar guiderate_fraction = local_fraction_lambda(guided_group); + // we add a factor here to avoid switching due to numerical instability + const Scalar factor = 1.01; + if (current_rate_fraction > (guiderate_fraction * factor)) { + return std::make_pair(true, guiderate_fraction / current_rate_fraction); + } + } + } + } + + if (this->schedule_.hasWell(name) && this->wellState().well(name).group_target) { + // for wells we already have computed the target + + // Switch sign since 'rates' are negative for producers. + const Scalar current_well_rate_available = -tcalc.calcModeRateFromRates(rates); + const Scalar group_target_rate_available = *this->wellState().well(name).group_target; + Scalar scale = 1.0; + if (current_well_rate_available > 1e-12) { + scale = group_target_rate_available / current_well_rate_available; + } + + return std::make_pair(current_well_rate_available > group_target_rate_available, scale); + } + + // we need to find out the level where the current well is applied to the local reduction + std::size_t local_reduction_level = 0; + for (std::size_t ii = 1; ii < num_ancestors; ++ii) { + const int num_gr_ctrl = this->groupControlledWells(chain[ii], + /*always_included_child=*/"", + /*is_producer=*/true, + /*injection_phase_not_used=*/Phase::OIL); + if (this->guide_rate_.has(chain[ii]) && num_gr_ctrl > 0) { + local_reduction_level = ii; + } + } + + // Compute portion of target corresponding to current_rate_available + Scalar target = orig_target; + for (std::size_t ii = 0; ii < num_ancestors; ++ii) { + if ((ii == 0) || this->guide_rate_.has(chain[ii])) { + // Apply local reductions only at the control level + // (top) and for levels where we have a specified + // group guide rate. + if (local_reduction_level >= ii) { + target -= local_reduction_lambda(chain[ii]); + } + // Add my reduction back at the level where it is included in the local reduction + if (local_reduction_level == ii) { + target += current_rate_available * efficiency_factor; + } + } + target *= local_fraction_lambda(chain[ii + 1]); + } + // Avoid negative target rates comming from too large local reductions. + const Scalar target_rate_available = std::max(Scalar(1e-12), target / efficiency_factor); + + Scalar scale = 1.0; + if (current_rate_available > 1e-12) + scale = target_rate_available / current_rate_available; + + return std::make_pair(current_rate_available > target_rate_available, scale); +} + +template +std::map +WellGroupHelper::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 +WellGroupHelper::getGuideRate(const std::string& name, + const GuideRateModel::Target target) const +{ + if (this->schedule_.hasWell(name, this->report_step_)) { + if (this->guide_rate_.has(name) || this->guide_rate_.hasPotentials(name)) { + return this->guide_rate_.get(name, target, this->getWellRateVector(name)); + } else { + return 0.0; + } + } + + if (this->guide_rate_.has(name)) { + return this->guide_rate_.get(name, target, this->getProductionGroupRateVector(name)); + } + + Scalar total_guide_rate = 0.0; + const Group& group = this->schedule_.getGroup(name, this->report_step_); + + for (const std::string& group_name : group.groups()) { + const Group::ProductionCMode& current_group_control + = this->groupState().production_control(group_name); + if (current_group_control == Group::ProductionCMode::FLD + || current_group_control == Group::ProductionCMode::NONE) { + // accumulate from sub wells/groups + total_guide_rate += this->getGuideRate(group_name, target); + } + } + + for (const std::string& well_name : group.wells()) { + const auto& well_tmp = this->schedule_.getWell(well_name, this->report_step_); + + if (well_tmp.isInjector()) + continue; + + const auto well_index = this->wellState().index(well_name); + if (!well_index.has_value()) + continue; + + if (!this->wellState().wellIsOwned(well_index.value(), well_name)) // Only sum once + { + continue; + } + + const auto& ws = this->wellState().well(well_index.value()); + if (ws.status == Well::Status::SHUT) + continue; + + if (!this->wellState().isProductionGrup(well_name)) + continue; + + // Only count wells under group control or the ru + if (!this->wellState().isProductionGrup(well_name)) + continue; + + total_guide_rate += this->getGuideRate(well_name, target); + } + return total_guide_rate; +} + +template +GuideRate::RateVector +WellGroupHelper::getProductionGroupRateVector(const std::string& group_name) const +{ + return this->getGuideRateVector_(this->groupState().production_rates(group_name)); +} + +template +Scalar +WellGroupHelper::getWellGroupTargetInjector(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + Phase injection_phase, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + DeferredLogger& deferred_logger) const +{ + // This function computes a wells group target. + // 'parent' will be the name of 'group'. But if we recurse, 'name' and + // 'parent' will stay fixed while 'group' will be higher up + // in the group tree. + // Eficiencyfactor is the well efficiency factor for the first group the well is + // part of. Later it is the accumulated factor including the group efficiency factor + // of the child of group. + auto current_group_control = this->groupState().injection_control(group.name(), injection_phase); + if (current_group_control == Group::InjectionCMode::FLD + || current_group_control == Group::InjectionCMode::NONE) { + // Return if we are not available for parent group. + if (!group.injectionGroupControlAvailable(injection_phase)) { + return std::numeric_limits::max(); + } + // Otherwise: check production share of parent's control. + const auto& parent_group = this->schedule_.getGroup(group.parent(), this->report_step_); + return this->getWellGroupTargetInjector(name, + parent, + parent_group, + rates, + injection_phase, + efficiency_factor * group.getGroupEfficiencyFactor(), + resv_coeff, + deferred_logger); + } + + // This can be false for FLD-controlled groups, we must therefore + // check for FLD first (done above). + if (!group.isInjectionGroup()) { + return std::numeric_limits::max(); + } + + // If we are here, we are at the topmost group to be visited in the recursion. + // This is the group containing the control we will check against. + Scalar sales_target = 0; + if (this->schedule_[this->report_step_].gconsale().has(group.name())) { + const auto& gconsale + = this->schedule_[this->report_step_].gconsale().get(group.name(), this->summary_state_); + sales_target = gconsale.sales_target; + } + WGHelpers::InjectionTargetCalculator tcalc { + current_group_control, + this->phase_usage_info_, + resv_coeff, + group.name(), + sales_target, + this->groupState(), + injection_phase, + group.has_gpmaint_control(injection_phase, current_group_control), + deferred_logger}; + + WGHelpers::FractionCalculator fcalc {this->schedule_, + *this, + this->summary_state_, + this->report_step_, + &this->guide_rate_, + tcalc.guideTargetMode(), + /*is_production_group=*/false, + injection_phase}; + + auto local_fraction_lambda = [&](const std::string& child, const std::string& always_incluced_name) { + return fcalc.localFraction(child, always_incluced_name); + }; + + auto local_reduction_lambda = [&](const std::string& group_name) { + const std::vector& group_target_reductions + = this->groupState().injection_reduction_rates(group_name); + return tcalc.calcModeRateFromRates(group_target_reductions); + }; + + std::optional ctrl; + if (!group.has_gpmaint_control(injection_phase, current_group_control)) + ctrl = group.injectionControls(injection_phase, this->summary_state_); + + + const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); + // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. + // Then ... + // TODO finish explanation. + const Scalar current_rate_available + = tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. + const auto chain = this->groupChainTopBot(name, group.name()); + // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. + const std::size_t num_ancestors = chain.size() - 1; + // we need to find out the level where the local reduction is applied + std::size_t local_reduction_level = 0; + for (std::size_t ii = 1; ii < num_ancestors; ++ii) { + const int num_gr_ctrl = this->groupControlledWells(chain[ii], + /*always_included_child=*/"", + /*is_production_group=*/false, + injection_phase); + if (this->guide_rate_.has(chain[ii], injection_phase) && num_gr_ctrl > 0) { + local_reduction_level = ii; + } + } + + // Compute portion of target corresponding to current_rate_available + Scalar target = orig_target; + for (std::size_t ii = 0; ii < num_ancestors; ++ii) { + if ((ii == 0) || this->guide_rate_.has(chain[ii], injection_phase)) { + // Apply local reductions only at the control level + // (top) and for levels where we have a specified + // group guide rate. + if (local_reduction_level >= ii) { + target -= local_reduction_lambda(chain[ii]); + } + + // If we are under individual control we need to add the wells rate back at the level where it is + // included in the local reduction + if (local_reduction_level == ii && !this->wellState().isInjectionGrup(name)) { + target += current_rate_available * efficiency_factor; + } + } + if (!this->wellState().isInjectionGrup(name)) { + target *= local_fraction_lambda(chain[ii + 1], name); + } else { + target *= local_fraction_lambda(chain[ii + 1], chain[ii + 1]); + } + } + // Avoid negative target rates comming from too large local reductions. + return std::max(Scalar(0.0), target / efficiency_factor); +} + +template +Scalar +WellGroupHelper::getWellGroupTargetProducer(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + DeferredLogger& deferred_logger) const +{ + // This function computes a wells group target. + // 'parent' will be the name of 'group'. But if we recurse, 'name' and + // 'parent' will stay fixed while 'group' will be higher up + // in the group tree. + // Eficiencyfactor is the well efficiency factor for the first group the well is + // part of. Later it is the accumulated factor including the group efficiency factor + // of the child of group. + OPM_TIMEFUNCTION(); + const Group::ProductionCMode& current_group_control = this->groupState().production_control(group.name()); + + if (current_group_control == Group::ProductionCMode::FLD + || current_group_control == Group::ProductionCMode::NONE) { + // Return if we are not available for parent group. + if (!group.productionGroupControlAvailable()) { + return std::numeric_limits::max(); + } + // Otherwise: check production share of parent's control. + const auto& parent_group = this->schedule_.getGroup(group.parent(), this->report_step_); + return this->getWellGroupTargetProducer(name, + parent, + parent_group, + rates, + efficiency_factor * group.getGroupEfficiencyFactor(), + resv_coeff, + deferred_logger); + } + + // This can be false for FLD-controlled groups, we must therefore + // check for FLD first (done above). + if (!group.isProductionGroup()) { + return std::numeric_limits::max(); + } + + // If we are here, we are at the topmost group to be visited in the recursion. + // This is the group containing the control we will check against. + + // gconsale may adjust the grat target. + // the adjusted rates is send to the targetCalculator + Scalar grat_target_from_sales = 0.0; + if (this->groupState().has_grat_sales_target(group.name())) + grat_target_from_sales = this->groupState().grat_sales_target(group.name()); + + WGHelpers::TargetCalculator tcalc {current_group_control, + this->phase_usage_info_, + resv_coeff, + grat_target_from_sales, + group.name(), + this->groupState(), + group.has_gpmaint_control(current_group_control)}; + + WGHelpers::FractionCalculator fcalc {this->schedule_, + *this, + this->summary_state_, + this->report_step_, + &this->guide_rate_, + tcalc.guideTargetMode(), + true, + Phase::OIL}; + auto local_fraction_lambda = [&](const std::string& child, const std::string& always_incluced_name) { + return fcalc.localFraction(child, always_incluced_name); + }; + + auto local_reduction_lambda = [&](const std::string& group_name) { + const std::vector& group_target_reductions + = this->groupState().production_reduction_rates(group_name); + return tcalc.calcModeRateFromRates(group_target_reductions); + }; + + std::optional ctrl; + if (!group.has_gpmaint_control(current_group_control)) + ctrl = group.productionControls(this->summary_state_); + + const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); + // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. + // Then ... + // TODO finish explanation. + const Scalar current_rate_available + = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. + const auto chain = this->groupChainTopBot(name, group.name()); + // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. + const std::size_t num_ancestors = chain.size() - 1; + // we need to find out the level where the local reduction is applied + std::size_t local_reduction_level = 0; + for (std::size_t ii = 1; ii < num_ancestors; ++ii) { + const int num_gr_ctrl = this->groupControlledWells(chain[ii], + /*always_included_child=*/"", + /*is_production_group=*/true, + /*injection_phase=*/Phase::OIL); + if (this->guide_rate_.has(chain[ii]) && num_gr_ctrl > 0) { + local_reduction_level = ii; + } + } + // Compute portion of target corresponding to current_rate_available + Scalar target = orig_target; + for (std::size_t ii = 0; ii < num_ancestors; ++ii) { + if ((ii == 0) || this->guide_rate_.has(chain[ii])) { + // Apply local reductions only at the control level + // (top) and for levels where we have a specified + // group guide rate. + if (local_reduction_level >= ii) { + target -= local_reduction_lambda(chain[ii]); + } + // If we are under individual control we need to add the wells rate back at the level where it is + // included in the local reduction + if (local_reduction_level == ii && !this->wellState().isProductionGrup(name)) { + target += current_rate_available * efficiency_factor; + } + } + if (this->wellState().isProductionGrup(name)) { + target *= local_fraction_lambda(chain[ii + 1], chain[ii + 1]); + } else { + target *= local_fraction_lambda(chain[ii + 1], name); + } + } + // Avoid negative target rates comming from too large local reductions. + return std::max(Scalar(0.0), target / efficiency_factor); + ; +} + +template +GuideRate::RateVector +WellGroupHelper::getWellRateVector(const std::string& name) const +{ + return this->getGuideRateVector_(this->wellState().currentWellRates(name)); +} + +template +std::vector +WellGroupHelper::groupChainTopBot(const std::string& bottom, + const std::string& top) const +{ + // Get initial parent, 'bottom' can be a well or a group. + std::string parent; + if (this->schedule_.hasWell(bottom, this->report_step_)) { + parent = this->schedule_.getWell(bottom, this->report_step_).groupName(); + } else { + parent = this->schedule_.getGroup(bottom, this->report_step_).parent(); + } + + // Build the chain from bottom to top. + std::vector chain; + chain.push_back(bottom); + chain.push_back(parent); + while (parent != top) { + parent = this->schedule_.getGroup(parent, this->report_step_).parent(); + chain.push_back(parent); + } + assert(chain.back() == top); + + // Reverse order and return. + std::reverse(chain.begin(), chain.end()); + return chain; +} + +template +int +WellGroupHelper::groupControlledWells(const std::string& group_name, + const std::string& always_included_child, + const bool is_production_group, + const Phase injection_phase) const +{ + auto num_wells = is_production_group + ? this->groupState().number_of_wells_under_group_control(group_name) + : this->groupState().number_of_wells_under_inj_group_control(group_name, injection_phase); + if (this->schedule_.hasWell(always_included_child, this->report_step_)) { + const bool isInGroup = this->isInGroupChainTopBot_(always_included_child, group_name); + const bool already_included = is_production_group + ? this->wellState().isProductionGrup(always_included_child) + : this->wellState().isInjectionGrup(always_included_child); + if (!already_included && isInGroup) { + num_wells++; + } + } + return num_wells; +} + +template +void +WellGroupHelper::setCmodeGroup(const Group& group, GroupState& group_state) const +{ + + for (const std::string& group_name : group.groups()) { + this->setCmodeGroup(this->schedule_.getGroup(group_name, this->report_step_), group_state); + } + + // use NONE as default control + const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; + for (Phase phase : all) { + if (!group_state.has_injection_control(group.name(), phase)) { + group_state.injection_control(group.name(), phase, Group::InjectionCMode::NONE); + } + } + if (!group_state.has_production_control(group.name())) { + group_state.production_control(group.name(), Group::ProductionCMode::NONE); + } + + const auto& events = this->schedule_[this->report_step_].wellgroup_events(); + if (group.isInjectionGroup() && events.hasEvent(group.name(), ScheduleEvents::GROUP_INJECTION_UPDATE)) { + + for (Phase phase : all) { + if (!group.hasInjectionControl(phase)) + continue; + + const auto& controls = group.injectionControls(phase, this->summary_state_); + group_state.injection_control(group.name(), phase, controls.cmode); + } + } + + if (group.isProductionGroup() && events.hasEvent(group.name(), ScheduleEvents::GROUP_PRODUCTION_UPDATE)) { + const auto controls = group.productionControls(this->summary_state_); + group_state.production_control(group.name(), controls.cmode); + } + + if (group.has_gpmaint_control(Group::ProductionCMode::RESV)) { + group_state.production_control(group.name(), Group::ProductionCMode::RESV); + } + for (Phase phase : all) { + if (group.has_gpmaint_control(phase, Group::InjectionCMode::RATE)) { + group_state.injection_control(group.name(), phase, Group::InjectionCMode::RATE); + } else if (group.has_gpmaint_control(phase, Group::InjectionCMode::RESV)) { + group_state.injection_control(group.name(), phase, Group::InjectionCMode::RESV); + } + } + + if (this->schedule_[this->report_step_].gconsale().has(group.name())) { + group_state.injection_control(group.name(), Phase::GAS, Group::InjectionCMode::SALE); + } +} + +// NOTE: setRegionAveragePressureCalculator() is a template member function and must be +// defined in the header file (see WellGroupHelper.hpp). +// +// See: updateGpMaintTargetForGroups() for detailed rationale. + +template +Scalar +WellGroupHelper::sumSolventRates(const Group& group, const bool is_injector) const +{ + Scalar rate = 0.0; + for (const std::string& group_name : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + const auto& gefac = group_tmp.getGroupEfficiencyFactor(); + rate += gefac * this->sumSolventRates(group_tmp, is_injector); + } + + for (const std::string& well_name : group.wells()) { + const auto well_index = this->wellState().index(well_name); + if (!well_index.has_value()) + continue; + + if (!this->wellState().wellIsOwned(well_index.value(), well_name)) // Only sum once + { + continue; + } + + const auto& well_ecl = this->schedule_.getWell(well_name, this->report_step_); + // only count producers or injectors + if ((well_ecl.isProducer() && is_injector) || (well_ecl.isInjector() && !is_injector)) + continue; + + const auto& ws = this->wellState().well(well_index.value()); + + if (ws.status == Well::Status::SHUT) + continue; + + const Scalar factor + = well_ecl.getEfficiencyFactor() * this->wellState()[well_ecl.name()].efficiency_scaling_factor; + if (is_injector) + rate += factor * ws.sum_solvent_rates(); + else + rate -= factor * ws.sum_solvent_rates(); + } + return rate; +} + +template +Scalar +WellGroupHelper::sumWellPhaseRates(bool res_rates, + const Opm::Group& group, + const int phase_pos, + const bool injector, + const bool network) const +{ + // Only obtain satellite rates once (on rank 0) + Scalar rate = 0.0; + if (this->wellState().isRank0() && (group.hasSatelliteProduction() || group.hasSatelliteInjection())) { + if (injector) { + rate = this->satelliteInjectionRate_( + this->schedule_[this->report_step_], group, phase_pos, res_rates); + } else { + const auto rate_comp = this->selectRateComponent_(phase_pos); + if (rate_comp.has_value()) { + rate = this->satelliteProductionRate_( + this->schedule_[this->report_step_], group, *rate_comp, res_rates); + } + } + // Satellite groups have no sub groups/wells so we're done + return rate; + } + + for (const std::string& group_name : group.groups()) { + const auto& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + const auto& gefac = group_tmp.getGroupEfficiencyFactor(network); + rate += gefac * this->sumWellPhaseRates(res_rates, group_tmp, phase_pos, injector, network); + } + + for (const std::string& well_name : group.wells()) { + const auto well_index = this->wellState().index(well_name); + if (!well_index.has_value()) + continue; + + if (!this->wellState().wellIsOwned(well_index.value(), well_name)) // Only sum once + { + continue; + } + + const auto& well_ecl = this->schedule_.getWell(well_name, this->report_step_); + // only count producers or injectors + if ((well_ecl.isProducer() && injector) || (well_ecl.isInjector() && !injector)) + continue; + + const auto& ws = this->wellState().well(well_index.value()); + if (ws.status == Opm::Well::Status::SHUT) + continue; + + const Scalar factor = well_ecl.getEfficiencyFactor(network) + * this->wellState().well(well_index.value()).efficiency_scaling_factor; + if (res_rates) { + const auto& well_rates = ws.reservoir_rates; + if (injector) + rate += factor * well_rates[phase_pos]; + else + rate -= factor * well_rates[phase_pos]; + } else { + const auto& well_rates = ws.surface_rates; + if (injector) + rate += factor * well_rates[phase_pos]; + else + rate -= factor * well_rates[phase_pos]; + } + } + return rate; +} + +template +Scalar +WellGroupHelper::sumWellResRates(const Group& group, + const int phasePos, + const bool injector) const +{ + return this->sumWellPhaseRates(/*res_rates=*/true, group, phasePos, injector); +} + +template +Scalar +WellGroupHelper::sumWellSurfaceRates(const Group& group, + const int phase_pos, + const bool injector) const +{ + return this->sumWellPhaseRates(/*res_rates=*/false, group, phase_pos, injector); +} + +// NOTE: updateGpMaintTargetForGroups() is a template member function and must be +// defined in the header file (see WellGroupHelper.hpp). +// +// Rationale: +// - The RegionalValues template parameter depends on BlackoilWellModel's FluidSystem, +// which varies by derived class and is not available in this compilation unit. +// - Template functions must be visible at their point of instantiation, so the +// implementation cannot be in this .cpp file. +// - We cannot use explicit template instantiation because the specific RegionalValues +// type depends on types from derived classes (e.g., AverageRegionalPressure). +// +// Historical Note: +// - The original implementation in WellGroupHelpers.cpp (static utility class) could keep +// the implementation in the .cpp file because it hard-coded BlackOilFluidSystem +// and used explicit template instantiation for that specific type. +// - As a member-based class, WellGroupHelper must support whatever FluidSystem the TypeTag +// provides, so we cannot hard-code the type and must use header-based implementation. +// +// If moved here, the linker will report: +// undefined reference to +// `WellGroupHelper<...>::updateGpMaintTargetForGroups>(const Group&, ...)` + +template +int +WellGroupHelper::updateGroupControlledWells(const bool is_production_group, + const Phase injection_phase, + GroupState& group_state, + DeferredLogger& deferred_logger) const +{ + OPM_TIMEFUNCTION(); + const auto& group_name = "FIELD"; + return this->updateGroupControlledWellsRecursive_( + group_name, is_production_group, injection_phase, group_state, deferred_logger); +} + +template +void +WellGroupHelper::updateGroupProductionRates(const Group& group, + GroupState& group_state) const +{ + OPM_TIMEFUNCTION(); + for (const std::string& group_name : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + this->updateGroupProductionRates(group_tmp, group_state); + } + const int np = this->wellState().numPhases(); + std::vector rates(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + rates[phase] = this->sumWellPhaseRates( + /*res_rates=*/false, group, phase, /*injector=*/false, /*network=*/false); + } + group_state.update_production_rates(group.name(), rates); +} + +template +void +WellGroupHelper::updateGroupTargetReduction(const Group& group, + const bool is_injector, + GroupState& group_state) const +{ + OPM_TIMEFUNCTION(); + std::vector group_target_reduction(this->wellState().numPhases(), 0.0); + this->updateGroupTargetReductionRecursive_(group, is_injector, group_target_reduction, group_state); +} + +template +void +WellGroupHelper::updateNetworkLeafNodeProductionRates( + GroupState& group_state) const +{ + const auto& network = this->schedule_[this->report_step_].network(); + if (network.active()) { + const int np = this->wellState().numPhases(); + for (const auto& group_name : network.leaf_nodes()) { + std::vector network_rates(np, 0.0); + if (this->schedule_[this->report_step_].groups.has( + group_name)) { // Allow empty leaf nodes that are not groups + const auto& group = this->schedule_[this->report_step_].groups.get(group_name); + if (group.numWells() > 0) { + for (int phase = 0; phase < np; ++phase) { + network_rates[phase] = this->sumWellPhaseRates( + /*res_rates=*/false, group, phase, /*injector=*/false, /*network=*/true); + } + } + } + group_state.update_network_leaf_node_production_rates(group_name, network_rates); + } + } +} + +template +void +WellGroupHelper::updateREINForGroups(const Group& group, + bool sum_rank, + GroupState& group_state) const +{ + OPM_TIMEFUNCTION(); + const int np = this->wellState().numPhases(); + for (const std::string& groupName : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(groupName, this->report_step_); + this->updateREINForGroups(group_tmp, sum_rank, group_state); + } + + std::vector rein(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + rein[phase] = this->sumWellPhaseRates( + /*res_rates=*/false, group, phase, /*injector=*/false, /*network=*/false); + } + + // add import rate and subtract consumption rate for group for gas + if (sum_rank) { + const auto& pu = this->wellState().phaseUsageInfo(); + if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { + const auto& [consumption_rate, import_rate] = group_state.gconsump_rates(group.name()); + const auto ig = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); + rein[ig] += import_rate; + rein[ig] -= consumption_rate; + } + } + + group_state.update_injection_rein_rates(group.name(), rein); +} + +template +void +WellGroupHelper::updateReservoirRatesInjectionGroups( + const Group& group, GroupState& group_state) const +{ + OPM_TIMEFUNCTION(); + for (const std::string& group_name : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + this->updateReservoirRatesInjectionGroups(group_tmp, group_state); + } + const int np = this->wellState().numPhases(); + std::vector resv(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + resv[phase] = this->sumWellPhaseRates( + /*res_rates=*/true, group, phase, /*injector=*/true, /*network=*/false); + } + group_state.update_injection_reservoir_rates(group.name(), resv); +} + +template +void +WellGroupHelper::updateState(WellState& well_state, + GroupState& group_state) +{ + this->well_state_ = &well_state; + this->group_state_ = &group_state; +} + +template +void +WellGroupHelper::updateSurfaceRatesInjectionGroups(const Group& group, + GroupState& group_state) const +{ + OPM_TIMEFUNCTION(); + for (const std::string& group_name : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + this->updateSurfaceRatesInjectionGroups(group_tmp, group_state); + } + const int np = this->wellState().numPhases(); + std::vector rates(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + rates[phase] = this->sumWellPhaseRates( + /*res_rates=*/false, + group, + phase, + /*injector=*/true, + /*network=*/false); + } + group_state.update_injection_surface_rates(group.name(), rates); +} + +template +void +WellGroupHelper::updateVREPForGroups(const Group& group, + GroupState& group_state) const +{ + OPM_TIMEFUNCTION(); + for (const std::string& group_name : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + this->updateVREPForGroups(group_tmp, group_state); + } + const int np = this->wellState().numPhases(); + Scalar resv = 0.0; + for (int phase = 0; phase < np; ++phase) { + resv += this->sumWellPhaseRates( + /*res_rates=*/true, group, phase, /*injector=*/false, /*network=*/false); + } + group_state.update_injection_vrep_rate(group.name(), resv); +} + +template +void +WellGroupHelper::updateWellRates(const Group& group, + const WellState& well_state_nupcol, + WellState& well_state) const +{ + OPM_TIMEFUNCTION(); + for (const std::string& group_name : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + this->updateWellRates(group_tmp, well_state_nupcol, well_state); + } + const int np = well_state.numPhases(); + for (const std::string& well_name : group.wells()) { + std::vector rates(np, 0.0); + const auto well_index = well_state.index(well_name); + if (well_index.has_value()) { // the well is found on this node + const auto& well_tmp = this->schedule_.getWell(well_name, this->report_step_); + int sign = 1; + // production wellRates are negative. The users of currentWellRates uses the convention in + // opm-common that production and injection rates are positive. + if (!well_tmp.isInjector()) + sign = -1; + const auto& ws = well_state_nupcol.well(well_index.value()); + for (int phase = 0; phase < np; ++phase) { + rates[phase] = sign * ws.surface_rates[phase]; + } + } + well_state.setCurrentWellRates(well_name, rates); + } +} + +template +void +WellGroupHelper::updateWellRatesFromGroupTargetScale( + const Scalar scale, + const Group& group, + bool is_injector, + WellState& well_state) const +{ + OPM_TIMEFUNCTION(); + for (const std::string& group_name : group.groups()) { + bool individual_control = false; + if (is_injector) { + const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; + for (Phase phase : all) { + const Group::InjectionCMode& current_group_control + = this->groupState().injection_control(group_name, phase); + individual_control = individual_control + || (current_group_control != Group::InjectionCMode::FLD + && current_group_control != Group::InjectionCMode::NONE); + } + } else { + const Group::ProductionCMode& current_group_control + = this->groupState().production_control(group_name); + individual_control = (current_group_control != Group::ProductionCMode::FLD + && current_group_control != Group::ProductionCMode::NONE); + } + if (!individual_control) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + this->updateWellRatesFromGroupTargetScale(scale, group_tmp, is_injector, well_state); + } + } + + const int np = well_state.numPhases(); + for (const std::string& well_name : group.wells()) { + const auto& well_tmp = this->schedule_.getWell(well_name, this->report_step_); + + if (well_tmp.isProducer() && is_injector) + continue; + + if (well_tmp.isInjector() && !is_injector) + continue; + + const auto well_index = well_state.index(well_name); + if (!well_index.has_value()) + continue; + + auto& ws = well_state.well(well_index.value()); + if (ws.status == Well::Status::SHUT) + continue; + + // scale rates + if (is_injector) { + if (ws.injection_cmode == Well::InjectorCMode::GRUP) + for (int phase = 0; phase < np; phase++) { + ws.surface_rates[phase] *= scale; + } + } else { + if (ws.production_cmode == Well::ProducerCMode::GRUP) { + for (int phase = 0; phase < np; phase++) { + ws.surface_rates[phase] *= scale; + } + } + } + } +} + +template +std::pair, Scalar> +WellGroupHelper::worstOffendingWell(const Group& group, + const Group::ProductionCMode& offended_control, + const Parallel::Communication& comm, + DeferredLogger& deferred_logger) const +{ + std::pair, Scalar> offending_well {std::nullopt, 0.0}; + for (const std::string& child_group : group.groups()) { + const auto& this_group = this->schedule_.getGroup(child_group, this->report_step_); + const auto& offending_well_this = this->worstOffendingWell( + this_group, offended_control, comm, deferred_logger + ); + if (offending_well_this.second > offending_well.second) { + offending_well = offending_well_this; + } + } + + for (const std::string& child_well : group.wells()) { + + const auto& well_index = this->wellState().index(child_well); + Scalar violating_rate = 0.0; + Scalar prefered_rate = 0.0; + const auto& pu = this->phase_usage_info_; + if (well_index.has_value() && this->wellState().wellIsOwned(well_index.value(), child_well)) { + const auto& ws = this->wellState().well(child_well); + switch (offended_control) { + case Group::ProductionCMode::ORAT: { + const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + violating_rate = ws.surface_rates[oil_pos]; + break; + } + case Group::ProductionCMode::GRAT: { + const int gas_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); + violating_rate = ws.surface_rates[gas_pos]; + break; + } + case Group::ProductionCMode::WRAT: { + const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); + violating_rate = ws.surface_rates[water_pos]; + break; + } + case Group::ProductionCMode::LRAT: { + assert(pu.phaseIsActive(IndexTraits::oilPhaseIdx)); + assert(pu.phaseIsActive(IndexTraits::waterPhaseIdx)); + const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); + violating_rate = ws.surface_rates[oil_pos] + ws.surface_rates[water_pos]; + break; + } + case Group::ProductionCMode::RESV: + for (int p = 0; p < this->wellState().numPhases(); ++p) { + violating_rate += ws.reservoir_rates[p]; + } + break; + case Group::ProductionCMode::NONE: + break; + case Group::ProductionCMode::FLD: + break; + case Group::ProductionCMode::PRBL: + OPM_DEFLOG_THROW(std::runtime_error, + "Group " + group.name() + " GroupProductionCMode PRBL not implemented", + deferred_logger); + break; + case Group::ProductionCMode::CRAT: + OPM_DEFLOG_THROW(std::runtime_error, + "Group " + group.name() + " GroupProductionCMode CRAT not implemented", + deferred_logger); + break; + } + const auto preferred_phase + = this->schedule_.getWell(child_well, this->report_step_).getPreferredPhase(); + switch (preferred_phase) { + case Phase::OIL: { + const int oil_pos + = this->phase_usage_info_.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + prefered_rate = ws.surface_rates[oil_pos]; + break; + } + case Phase::GAS: { + const int gas_pos + = this->phase_usage_info_.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); + prefered_rate = ws.surface_rates[gas_pos]; + break; + } + case Phase::WATER: { + const int water_pos + = this->phase_usage_info_.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); + prefered_rate = ws.surface_rates[water_pos]; + break; + } + default: + // No others supported. + break; + } + } + violating_rate = comm.sum(violating_rate); + if (violating_rate < 0) { // only check producing wells + prefered_rate = comm.sum(prefered_rate); + Scalar fraction = prefered_rate < -1e-16 ? violating_rate / prefered_rate : 1.0; + if (fraction > offending_well.second) { + offending_well = {child_well, fraction}; + } + } + } + return offending_well; +} + +// --------------------------------------------------------------------- +// Private methods +// --------------------------------------------------------------------- + +template +std::string +WellGroupHelper::controlGroup_(const Group& group) const +{ + const Group::ProductionCMode& currentGroupControl = this->groupState().production_control(group.name()); + + if (currentGroupControl == Group::ProductionCMode::FLD + || currentGroupControl == Group::ProductionCMode::NONE) { + const auto& parent_name = group.control_group(); + if (parent_name) { + const auto& parent = this->schedule_.getGroup(parent_name.value(), this->report_step_); + return this->controlGroup_(parent); + } + } + + return group.name(); +} + +template +Opm::GuideRate::RateVector +WellGroupHelper::getGuideRateVector_(const std::vector& rates) const +{ + const auto& pu = this->phase_usage_info_; + Scalar oilRate = 0.0; + if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { + oilRate = rates[pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx)]; + } + + Scalar gasRate = 0.0; + if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { + gasRate = rates[pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx)]; + } + + Scalar waterRate = 0.0; + if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { + waterRate = rates[pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx)]; + } + + return {oilRate, gasRate, waterRate}; +} + +template +bool +WellGroupHelper::isInGroupChainTopBot_(const std::string& bottom, + const std::string& top) const +{ + // Get initial parent, 'bottom' can be a well or a group. + std::string parent; + if (this->schedule_.hasWell(bottom, this->report_step_)) { + parent = this->schedule_.getWell(bottom, this->report_step_).groupName(); + } else { + parent = this->schedule_.getGroup(bottom, this->report_step_).parent(); + } + + while (parent != top) { + parent = this->schedule_.getGroup(parent, this->report_step_).parent(); + if (parent == top) { + return true; + } else if (parent == "FIELD") { + return false; + } + } + return true; +} + +template +int +WellGroupHelper::phaseToActivePhaseIdx_(const Phase phase) const +{ + const auto& pu = this->phase_usage_info_; + int phase_pos = -1; + switch (phase) { + case Phase::GAS: + if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { + phase_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); + } + break; + case Phase::OIL: + if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { + phase_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + } + break; + case Phase::WATER: + if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { + phase_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); + } + break; + default: + // just to avoid warning + throw std::invalid_argument("unhandled phase enum"); + } + return phase_pos; +} + +template +Scalar +WellGroupHelper::satelliteInjectionRate_(const ScheduleState& sched, + const Group& group, + const int phase_pos, + bool res_rates) const +{ + Scalar rate = 0.0; + const auto& pu = this->phase_usage_info_; + if (group.hasSatelliteInjection()) { + std::optional ph; + for (const auto& [bo_phase, phase] : std::array {std::pair(IndexTraits::waterPhaseIdx, Phase::WATER), + std::pair(IndexTraits::oilPhaseIdx, Phase::OIL), + std::pair(IndexTraits::gasPhaseIdx, Phase::GAS)}) { + if (pu.phaseIsActive(bo_phase) && (pu.canonicalToActivePhaseIdx(bo_phase) == phase_pos)) { + ph = phase; + } + } + if (ph.has_value()) { + const auto& satellite_inj = sched.satelliteInjection(group.name()); + const auto& rate_ix = satellite_inj.rateIndex(ph.value()); + if (rate_ix.has_value()) { + const auto& satrates = satellite_inj[*rate_ix]; + if (!res_rates) { // surface rates + if (const auto& qs = satrates.surface(); qs.has_value()) { + rate = *qs; + } + // We don't support reservoir rates for satellite injection groups + } + } + } + } + return rate; +} + +template +Scalar +WellGroupHelper::satelliteProductionRate_( + const ScheduleState& sched, + const Group& group, + const GSatProd::GSatProdGroupProp::Rate rate_comp, + bool res_rates) const +{ + Scalar rate = 0.0; + if (group.hasSatelliteProduction()) { + const auto& gsat_prod = sched.gsatprod(); + if (!res_rates) { + rate = gsat_prod.get(group.name(), this->summary_state_).rate[rate_comp]; + } + // We don't support reservoir rates for satellite production groups + } + return rate; +} + +template +std::optional +WellGroupHelper::selectRateComponent_(const int phase_pos) const +{ + // TODO: this function can be wrong, phase_pos is not used anymore, this function requries checking and + // refactoring. + const auto& pu = this->phase_usage_info_; + using Rate = GSatProd::GSatProdGroupProp::Rate; + + for (const auto& [phase, rate_comp] : std::array {std::pair {IndexTraits::waterPhaseIdx, Rate::Water}, + std::pair {IndexTraits::oilPhaseIdx, Rate::Oil}, + std::pair {IndexTraits::gasPhaseIdx, Rate::Gas}}) { + if (pu.phaseIsActive(phase)) { + const auto active_phase_pos = pu.canonicalToActivePhaseIdx(phase); + if (active_phase_pos == phase_pos) { + return rate_comp; + } + } + } + + return std::nullopt; +} + +template +int +WellGroupHelper::updateGroupControlledWellsRecursive_( + const std::string& group_name, + const bool is_production_group, + const Phase injection_phase, + GroupState& group_state, + DeferredLogger& deferred_logger) const +{ + const Group& group = this->schedule_.getGroup(group_name, this->report_step_); + int num_wells = 0; + for (const std::string& child_group : group.groups()) { + + bool included = false; + if (is_production_group) { + const auto ctrl = group_state.production_control(child_group); + included + = included || (ctrl == Group::ProductionCMode::FLD) || (ctrl == Group::ProductionCMode::NONE); + } else { + const auto ctrl = group_state.injection_control(child_group, injection_phase); + included + = included || (ctrl == Group::InjectionCMode::FLD) || (ctrl == Group::InjectionCMode::NONE); + } + + if (included) { + num_wells += this->updateGroupControlledWellsRecursive_( + child_group, is_production_group, injection_phase, group_state, deferred_logger); + } else { + this->updateGroupControlledWellsRecursive_( + child_group, is_production_group, injection_phase, group_state, deferred_logger); + } + } + for (const std::string& child_well : group.wells()) { + bool included = false; + const Well& well = this->schedule_.getWell(child_well, this->report_step_); + if (is_production_group && well.isProducer()) { + included = included || this->wellState().isProductionGrup(child_well) || group.as_choke(); + } else if (!is_production_group && !well.isProducer()) { + const auto& well_controls = well.injectionControls(this->summary_state_); + auto injectorType = well_controls.injector_type; + if ((injection_phase == Phase::WATER && injectorType == InjectorType::WATER) + || (injection_phase == Phase::OIL && injectorType == InjectorType::OIL) + || (injection_phase == Phase::GAS && injectorType == InjectorType::GAS)) { + included = included || this->wellState().isInjectionGrup(child_well); + } + } + const auto ctrl1 = group_state.production_control(group.name()); + if (group.as_choke() + && ((ctrl1 == Group::ProductionCMode::FLD) || (ctrl1 == Group::ProductionCMode::NONE))) { + // The auto choke group has not own group control but inherits control from an ancestor group. + // Number of wells should be calculated as zero when wells of auto choke group do not deliver + // target. This behaviour is then similar to no-autochoke group with wells not on GRUP control. + // The rates of these wells are summed up. The parent group target is reduced with this rate. + // This reduced target becomes the target of the other child group of this parent. + const auto num_phases = this->wellState().numPhases(); + std::vector rates(num_phases, 0.0); + for (int phase_pos = 0; phase_pos < num_phases; ++phase_pos) { + rates[phase_pos] = this->sumWellSurfaceRates(group, phase_pos, /*injector=*/false); + } + + // Get the ancestor of the auto choke group that has group control (cmode != FLD, NONE) + const auto& control_group_name = this->controlGroup_(group); + const auto& control_group = this->schedule_.getGroup(control_group_name, this->report_step_); + const auto& ctrl = control_group.productionControls(this->summary_state_); + const auto& control_group_cmode = ctrl.cmode; + + const auto& group_guide_rate = group.productionControls(this->summary_state_).guide_rate; + + if (group_guide_rate > 0) { + // Guide rate is not default for the auto choke group + Scalar grat_target_from_sales = 0.0; + if (group_state.has_grat_sales_target(control_group_name)) + grat_target_from_sales = group_state.grat_sales_target(control_group_name); + + std::vector resv_coeff(num_phases, 1.0); + WGHelpers::TargetCalculator tcalc { + control_group_cmode, + this->phase_usage_info_, + resv_coeff, + grat_target_from_sales, + group.name(), + group_state, + group.has_gpmaint_control(control_group_cmode)}; + const auto& control_group_target = tcalc.groupTarget(ctrl, deferred_logger); + + // Calculates the guide rate of the parent group with control. + // It is allowed that the guide rate of this group is defaulted. The guide rate will be + // derived from the children groups + const auto& control_group_guide_rate + = this->getGuideRate(control_group_name, tcalc.guideTargetMode()); + + if (control_group_guide_rate > 0) { + // Target rate for the auto choke group + const Scalar target_rate + = control_group_target * group_guide_rate / control_group_guide_rate; + const Scalar current_rate = tcalc.calcModeRateFromRates(rates); + + if (current_rate < target_rate) + included = false; + } + } + } + + if (included) { + ++num_wells; + } + } + if (is_production_group) { + group_state.update_number_of_wells_under_group_control(group_name, num_wells); + } else { + group_state.update_number_of_wells_under_inj_group_control(group_name, injection_phase, num_wells); + } + + return num_wells; +} + +template +void +WellGroupHelper::updateGroupTargetReductionRecursive_( + const Group& group, + const bool is_injector, + std::vector& group_target_reduction, + GroupState& group_state) const +{ + const int np = this->wellState().numPhases(); + for (const std::string& sub_group_name : group.groups()) { + std::vector sub_group_target_reduction(np, 0.0); + const Group& sub_group = this->schedule_.getGroup(sub_group_name, this->report_step_); + this->updateGroupTargetReductionRecursive_( + sub_group, is_injector, sub_group_target_reduction, group_state); + + const Scalar sub_group_efficiency = sub_group.getGroupEfficiencyFactor(); + + // accumulate group contribution from sub group + if (is_injector) { + const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; + for (Phase phase : all) { + const auto phase_pos = this->phaseToActivePhaseIdx_(phase); + // the phase is not present + if (phase_pos == -1) + continue; + const Group::InjectionCMode& current_group_control + = group_state.injection_control(sub_group.name(), phase); + const bool individual_control = (current_group_control != Group::InjectionCMode::FLD + && current_group_control != Group::InjectionCMode::NONE); + const int num_group_controlled_wells + = this->groupControlledWells(sub_group.name(), + /*always_included_child=*/"", + !is_injector, + phase); + if (individual_control || num_group_controlled_wells == 0) { + group_target_reduction[phase_pos] += sub_group_efficiency + * this->sumWellSurfaceRates(sub_group, phase_pos, is_injector); + } else { + // Accumulate from this subgroup only if no group guide rate is set for it. + if (!this->guide_rate_.has(sub_group.name(), phase)) { + group_target_reduction[phase_pos] + += sub_group_efficiency * sub_group_target_reduction[phase_pos]; + } + } + } + } else { + const Group::ProductionCMode& current_group_control + = group_state.production_control(sub_group.name()); + const bool individual_control = (current_group_control != Group::ProductionCMode::FLD + && current_group_control != Group::ProductionCMode::NONE); + const int num_group_controlled_wells + = this->groupControlledWells(sub_group.name(), + /*always_included_child=*/"", + !is_injector, + /*injection_phase_not_used=*/Phase::OIL); + if (individual_control || num_group_controlled_wells == 0) { + for (int phase = 0; phase < np; phase++) { + group_target_reduction[phase] + += sub_group_efficiency * this->sumWellSurfaceRates(sub_group, phase, is_injector); + } + } else { + // The subgroup may participate in group control. + if (!this->guide_rate_.has(sub_group.name())) { + // Accumulate from this subgroup only if no group guide rate is set for it. + for (int phase = 0; phase < np; phase++) { + group_target_reduction[phase] + += sub_group_efficiency * sub_group_target_reduction[phase]; + } + } + } + } + } + + for (const std::string& well_name : group.wells()) { + const auto& well_tmp = this->schedule_.getWell(well_name, this->report_step_); + + if (well_tmp.isProducer() && is_injector) + continue; + + if (well_tmp.isInjector() && !is_injector) + continue; + + const auto well_index = this->wellState().index(well_name); + if (!well_index.has_value()) + continue; + + if (!this->wellState().wellIsOwned(well_index.value(), well_name)) // Only sum once + { + continue; + } + + const auto& ws = this->wellState().well(well_index.value()); + if (ws.status == Well::Status::SHUT) + continue; + + const Scalar efficiency = well_tmp.getEfficiencyFactor() + * this->wellState().well(well_index.value()).efficiency_scaling_factor; + + // add contribution from wells not under group control + if (is_injector) { + if (ws.injection_cmode != Well::InjectorCMode::GRUP) + for (int phase = 0; phase < np; phase++) { + group_target_reduction[phase] += ws.surface_rates[phase] * efficiency; + } + } else { + if ((ws.production_cmode != Well::ProducerCMode::GRUP)) { + if (!group.as_choke()) { + for (int phase = 0; phase < np; phase++) { + group_target_reduction[phase] -= ws.surface_rates[phase] * efficiency; + } + } + } + } + } + if (is_injector) { + group_state.update_injection_reduction_rates(group.name(), group_target_reduction); + } else { + group_state.update_production_reduction_rates(group.name(), group_target_reduction); + } +} + + +template class WellGroupHelper; +#ifdef FLOW_INSTANTIATE_FLOAT +template class WellGroupHelper; +#endif +} // namespace Opm diff --git a/opm/simulators/wells/WellGroupHelper.hpp b/opm/simulators/wells/WellGroupHelper.hpp new file mode 100644 index 0000000000..331403c447 --- /dev/null +++ b/opm/simulators/wells/WellGroupHelper.hpp @@ -0,0 +1,450 @@ +/* + Copyright 2025 Equinor ASA + + 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_WELLGROUPHELPER_HEADER_INCLUDED +#define OPM_WELLGROUPHELPER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +namespace Opm +{ + +template +class WellGroupHelper +{ +public: + // RAII guard for temporarily setting wellstate pointer + class WellStateGuard + { + public: + WellStateGuard(WellGroupHelper& wgHelper, WellState& well_state) + : wgHelper_ {wgHelper} + , previous_state_ptr_ {wgHelper_.well_state_} + { + // Set the new state directly + wgHelper_.well_state_ = &well_state; + } + + ~WellStateGuard() + { + // Restore the previous state + wgHelper_.well_state_ = previous_state_ptr_; + } + + // Delete copy and move operations + WellStateGuard(const WellStateGuard&) = delete; + WellStateGuard& operator=(const WellStateGuard&) = delete; + WellStateGuard(WellStateGuard&&) = delete; + WellStateGuard& operator=(WellStateGuard&&) = delete; + + private: + WellGroupHelper& wgHelper_; + const WellState* previous_state_ptr_; + }; + + // RAII guard for temporarily setting groupstate pointer + class GroupStateGuard + { + public: + GroupStateGuard(WellGroupHelper& wgHelper, GroupState& group_state) + : wgHelper_ {wgHelper} + , previous_state_ptr_ {wgHelper_.group_state_} + { + // Set the new state directly + wgHelper_.group_state_ = &group_state; + } + + ~GroupStateGuard() + { + // Restore the previous state + wgHelper_.group_state_ = previous_state_ptr_; + } + + // Delete copy and move operations + GroupStateGuard(const GroupStateGuard&) = delete; + GroupStateGuard& operator=(const GroupStateGuard&) = delete; + GroupStateGuard(GroupStateGuard&&) = delete; + GroupStateGuard& operator=(GroupStateGuard&&) = delete; + + private: + WellGroupHelper& wgHelper_; + const GroupState* previous_state_ptr_; + }; + + WellGroupHelper(WellState& well_state, + GroupState& group_state, + const Schedule& schedule, + const SummaryState& summary_state, + const GuideRate& guide_rate, + const PhaseUsageInfo& phase_usage_info); + + void accumulateGroupEfficiencyFactor(const Group& group, Scalar& factor) const; + + std::pair checkGroupConstraintsInj(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + const Phase injection_phase, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + const bool check_guide_rate, + DeferredLogger& deferred_logger) const; + + std::pair checkGroupConstraintsProd(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + 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; + + Scalar getWellGroupTargetInjector(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + const Phase injection_phase, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + DeferredLogger& deferred_logger) const; + + Scalar getWellGroupTargetProducer(const std::string& name, + const std::string& parent, + const Group& group, + const Scalar* rates, + const Scalar efficiency_factor, + const std::vector& resv_coeff, + DeferredLogger& deferred_logger) const; + + GuideRate::RateVector getWellRateVector(const std::string& name) const; + + std::vector groupChainTopBot(const std::string& bottom, const std::string& top) const; + + /// returns the number of wells that are actively under group control for a given group with name given + /// by group_name + int groupControlledWells(const std::string& group_name, + const std::string& always_included_child, + const bool is_production_group, + const Phase injection_phase) const; + + const GroupState& groupState() const + { + return *this->group_state_; + } + + const PhaseUsageInfo& phaseUsageInfo() const + { + return this->phase_usage_info_; + } + + GroupStateGuard pushGroupState(GroupState& group_state) + { + return GroupStateGuard(*this, group_state); + } + + WellStateGuard pushWellState(WellState& well_state) + { + return WellStateGuard(*this, well_state); + } + + void setCmodeGroup(const Group& group, GroupState& group_state) const; + + template + void + setRegionAveragePressureCalculator(const Group& group, + const FieldPropsManager& fp, + std::map>& + regional_average_pressure_calculator) const; + + void setReportStep(int report_step) + { + report_step_ = report_step; + } + + Scalar sumSolventRates(const Group& group, const bool is_injector) const; + + Scalar sumWellResRates(const Group& group, const int phase_pos, const bool injector) const; + + Scalar sumWellSurfaceRates(const Group& group, const int phase_pos, const bool injector) const; + + Scalar sumWellPhaseRates(bool res_rates, + const Group& group, + const int phase_pos, + const bool injector, + const bool network = false) const; + + template + void updateGpMaintTargetForGroups(const Group& group, + const RegionalValues& regional_values, + const double dt, + GroupState& group_state) const; + + /// update the number of wells that are actively under group control for a given group with name given by + /// group_name its main usage is to detect cases where there is no wells under group control + int updateGroupControlledWells(const bool is_production_group, + const Phase injection_phase, + GroupState& group_state, + DeferredLogger& deferred_logger) const; + + void updateGroupProductionRates(const Group& group, GroupState& group_state) const; + + void updateGroupTargetReduction(const Group& group, + const bool is_injector, + GroupState& group_state) const; + + void updateNetworkLeafNodeProductionRates(GroupState& group_state) const; + + void updateREINForGroups(const Group& group, bool sum_rank, GroupState& group_state) const; + + void updateReservoirRatesInjectionGroups(const Group& group, GroupState& group_state) const; + + void updateVREPForGroups(const Group& group, GroupState& group_state) const; + + void updateState(WellState& well_state, GroupState& group_state); + + void updateSurfaceRatesInjectionGroups(const Group& group, GroupState& group_state) const; + + void updateWellRates(const Group& group, + const WellState& well_state_nupcol, + WellState& well_state) const; + + const WellState& wellState() const + { + return *this->well_state_; + } + + void updateWellRatesFromGroupTargetScale(const Scalar scale, + const Group& group, + bool is_injector, + WellState& well_state) const; + + /// Returns the name of the worst offending well and its fraction (i.e. violated_phase / preferred_phase) + std::pair, Scalar> + worstOffendingWell(const Group& group, + const Group::ProductionCMode& offended_control, + const Parallel::Communication& comm, + DeferredLogger& deferred_logger) const; + +private: + std::string controlGroup_(const Group& group) const; + + GuideRate::RateVector getGuideRateVector_(const std::vector& rates) const; + + /// check if well/group bottom is a sub well/group of the group top + bool isInGroupChainTopBot_(const std::string& bottom, const std::string& top) const; + + int phaseToActivePhaseIdx_(const Phase phase) const; + + Scalar satelliteInjectionRate_(const ScheduleState& sched, + const Group& group, + const int phase_pos, + bool res_rates) const; + + Scalar satelliteProductionRate_(const ScheduleState& sched, + const Group& group, + const GSatProd::GSatProdGroupProp::Rate rate_comp, + bool res_rates) const; + + std::optional selectRateComponent_(const int phase_pos) const; + + int updateGroupControlledWellsRecursive_(const std::string& group_name, + const bool is_production_group, + const Phase injection_phase, + GroupState& group_state, + DeferredLogger& deferred_logger) const; + + void updateGroupTargetReductionRecursive_(const Group& group, + const bool is_injector, + std::vector& group_target_reduction, + GroupState& group_state) const; + + const WellState* well_state_ {nullptr}; + const GroupState* group_state_ {nullptr}; + const Schedule& schedule_; + const SummaryState& summary_state_; + const GuideRate& guide_rate_; + int report_step_ {0}; + // NOTE: The phase usage info seems to be read-only throughout the simulation, so it should be safe + // to store a reference to it here. + const PhaseUsageInfo& phase_usage_info_; +}; + +// ----------------------------------------------------------------------------- +// Template member function implementations +// ----------------------------------------------------------------------------- + +// NOTE: This template member function is defined here in the header because the +// AverageRegionalPressureType type parameter depends on derived class types that are +// not available when WellGroupHelper.cpp is compiled. Template functions +// must be visible at their instantiation point. +// See WellGroupHelper.cpp for detailed rationale. +template +template +void +WellGroupHelper::setRegionAveragePressureCalculator( + const Group& group, + const FieldPropsManager& fp, + std::map>& regional_average_pressure_calculator) + const +{ + for (const std::string& groupName : group.groups()) { + this->setRegionAveragePressureCalculator(this->schedule_.getGroup(groupName, this->report_step_), + fp, + regional_average_pressure_calculator); + } + const auto& gpm = group.gpmaint(); + if (!gpm) + return; + + const auto& reg = gpm->region(); + if (!reg) + return; + + if (regional_average_pressure_calculator.count(reg->first) == 0) { + const std::string name = (reg->first.rfind("FIP", 0) == 0) ? reg->first : "FIP" + reg->first; + const auto& fipnum = fp.get_int(name); + regional_average_pressure_calculator[reg->first] + = std::make_unique(fipnum); + } +} + +// NOTE: This template member function is defined in the header because the +// RegionalValues type parameter depends on derived class types that are +// not available when WellGroupHelper.cpp is compiled. Template functions +// must be visible at their instantiation point. +// See WellGroupHelper.cpp for detailed rationale. +template +template +void +WellGroupHelper::updateGpMaintTargetForGroups(const Group& group, + const RegionalValues& regional_values, + const double dt, + GroupState& group_state) const +{ + OPM_TIMEFUNCTION(); + for (const std::string& group_name : group.groups()) { + const Group& group_tmp = this->schedule_.getGroup(group_name, this->report_step_); + this->updateGpMaintTargetForGroups(group_tmp, regional_values, dt, group_state); + } + const auto& gpm = group.gpmaint(); + if (!gpm) + return; + + const auto& region = gpm->region(); + if (!region) + return; + + const auto [name, number] = *region; + const Scalar error = gpm->pressure_target() - regional_values.at(name)->pressure(number); + Scalar current_rate = 0.0; + const auto& pu = this->phase_usage_info_; + bool injection = true; + Scalar sign = 1.0; + switch (gpm->flow_target()) { + case GPMaint::FlowTarget::RESV_PROD: { + current_rate = -this->groupState().injection_vrep_rate(group.name()); + injection = false; + sign = -1.0; + break; + } + case GPMaint::FlowTarget::RESV_OINJ: { + if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { + const auto io = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + current_rate = group_state.injection_reservoir_rates(group.name())[io]; + } + break; + } + case GPMaint::FlowTarget::RESV_WINJ: { + if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { + const auto iw = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); + current_rate = group_state.injection_reservoir_rates(group.name())[iw]; + } + break; + } + case GPMaint::FlowTarget::RESV_GINJ: { + if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { + const auto ig = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); + current_rate = group_state.injection_reservoir_rates(group.name())[ig]; + } + break; + } + case GPMaint::FlowTarget::SURF_OINJ: { + if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { + const auto io = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); + current_rate = group_state.injection_surface_rates(group.name())[io]; + } + break; + } + case GPMaint::FlowTarget::SURF_WINJ: { + if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { + const auto iw = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); + current_rate = group_state.injection_surface_rates(group.name())[iw]; + } + break; + } + case GPMaint::FlowTarget::SURF_GINJ: { + if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { + const auto ig = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); + current_rate = group_state.injection_surface_rates(group.name())[ig]; + } + break; + } + default: + throw std::invalid_argument("Invalid Flow target type in GPMAINT"); + } + auto& gpmaint_state = group_state.gpmaint(group.name()); + // we only activate gpmaint if pressure is lower than the target regional pressure for injectors + // (i.e. error > 0) and higher for producers. + bool activate = (injection && error > 0) || (!injection && error < 0); + Scalar rate = 0.0; + if (activate) { + rate = gpm->rate(gpmaint_state, current_rate, error, dt); + } else { + gpm->resetState(gpmaint_state); + } + group_state.update_gpmaint_target(group.name(), std::max(Scalar {0.0}, sign * rate)); +} + +} // namespace Opm + +#endif diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp deleted file mode 100644 index 0279704534..0000000000 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ /dev/null @@ -1,2290 +0,0 @@ -/* - Copyright 2019 Norce. - Copyright 2020 Equinor ASA. - - 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 -#include -#include -#include - -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -namespace { - template - Opm::GuideRate::RateVector - getGuideRateVector(const std::vector& rates, const Opm::PhaseUsageInfo& pu) - { - Scalar oilRate = 0.0; - if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { - oilRate = rates[pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx)]; - } - - Scalar gasRate = 0.0; - if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { - gasRate = rates[pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx)]; - } - - Scalar waterRate = 0.0; - if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { - waterRate = rates[pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx)]; - } - - return {oilRate, gasRate, waterRate}; - } - -} // namespace Anonymous - -namespace Opm { - - template - Scalar WellGroupHelpers:: - sumWellPhaseRates(bool res_rates, - const Opm::Group& group, - const Opm::Schedule& schedule, - const Opm::WellState& wellState, - const SummaryState& summaryState, - const int reportStepIdx, - const int phasePos, - const bool injector, - const bool network) - { - // Only obtain satellite rates once (on rank 0) - Scalar rate = 0.0; - if (wellState.isRank0() && (group.hasSatelliteProduction() || group.hasSatelliteInjection())) { - if (injector) { - rate = satelliteInjectionRate(schedule[reportStepIdx], group, wellState.phaseUsageInfo(), phasePos, res_rates); - } else { - const auto rateComp = selectRateComponent(wellState.phaseUsageInfo(), phasePos); - if (rateComp.has_value()) { - rate = satelliteProductionRate(summaryState, schedule[reportStepIdx], group, *rateComp, res_rates); - } - } - // Satellite groups have no sub groups/wells so we're done - return rate; - } - - for (const std::string& groupName : group.groups()) { - const auto& groupTmp = schedule.getGroup(groupName, reportStepIdx); - const auto& gefac = groupTmp.getGroupEfficiencyFactor(network); - rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, summaryState, reportStepIdx, phasePos, injector, network); - } - - for (const std::string& wellName : group.wells()) { - const auto well_index = wellState.index(wellName); - if (!well_index.has_value()) - continue; - - if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once - { - continue; - } - - const auto& wellEcl = schedule.getWell(wellName, reportStepIdx); - // only count producers or injectors - if ((wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector)) - continue; - - const auto& ws = wellState.well(well_index.value()); - if (ws.status == Opm::Well::Status::SHUT) - continue; - - const Scalar factor = wellEcl.getEfficiencyFactor(network) * wellState[wellEcl.name()].efficiency_scaling_factor; - if (res_rates) { - const auto& well_rates = ws.reservoir_rates; - if (injector) - rate += factor * well_rates[phasePos]; - else - rate -= factor * well_rates[phasePos]; - } else { - const auto& well_rates = ws.surface_rates; - if (injector) - rate += factor * well_rates[phasePos]; - else - rate -= factor * well_rates[phasePos]; - } - } - return rate; - } - -template -Scalar WellGroupHelpers:: -satelliteInjectionRate(const ScheduleState& sched, - const Group& group, - const PhaseUsageInfo& pu, - const int phase_pos, - bool res_rates) -{ - Scalar rate = 0.0; - if (group.hasSatelliteInjection()) { - std::optional ph; - for (const auto& [bo_phase, phase] : std::array { - std::pair( IndexTraits::waterPhaseIdx, Phase::WATER), - std::pair( IndexTraits::oilPhaseIdx, Phase::OIL), - std::pair( IndexTraits::gasPhaseIdx, Phase::GAS) }) - { - if (pu.phaseIsActive(bo_phase) && (pu.canonicalToActivePhaseIdx(bo_phase) == phase_pos)) { - ph = phase; - } - } - if (ph.has_value()) { - const auto& satellite_inj = sched.satelliteInjection(group.name()); - const auto& rate_ix = satellite_inj.rateIndex(ph.value()); - if (rate_ix.has_value()) { - const auto& satrates = satellite_inj[*rate_ix]; - if (!res_rates) { // surface rates - if (const auto& qs = satrates.surface(); qs.has_value()) { - rate = *qs; - } - // We don't support reservoir rates for satellite injection groups - } - } - } - } - return rate; -} - -template -Scalar WellGroupHelpers:: -satelliteProductionRate(const SummaryState& summaryState, - const ScheduleState& sched, - const Group& group, - const GSatProd::GSatProdGroupProp::Rate rateComp, - bool res_rates) -{ - Scalar rate = 0.0; - if (group.hasSatelliteProduction()) { - const auto& gsatProd = sched.gsatprod(); - if (!res_rates) { - rate = gsatProd.get(group.name(), summaryState).rate[rateComp]; - } - // We don't support reservoir rates for satellite production groups - } - return rate; -} - -template -std::optional -WellGroupHelpers:: -selectRateComponent(const PhaseUsageInfo& pu, const int phasePos) -{ - // TODO: this function can be wrong, phasePos is not used anymore, this function requries checking and refactoring. - using Rate = GSatProd::GSatProdGroupProp::Rate; - - for (const auto& [phase, rateComp] : std::array { - std::pair { IndexTraits::waterPhaseIdx, Rate::Water }, - std::pair { IndexTraits::oilPhaseIdx, Rate::Oil }, - std::pair { IndexTraits::gasPhaseIdx, Rate::Gas } }) - { - if (pu.phaseIsActive(phase)) { - const auto activePhasePos = pu.canonicalToActivePhaseIdx(phase); - if (activePhasePos == phasePos) { - return rateComp; - } - } - } - - return std::nullopt; -} - -template -void WellGroupHelpers:: -setCmodeGroup(const Group& group, - const Schedule& schedule, - const SummaryState& summaryState, - const int reportStepIdx, - GroupState& group_state) -{ - - for (const std::string& groupName : group.groups()) { - setCmodeGroup(schedule.getGroup(groupName, reportStepIdx), - schedule, summaryState, reportStepIdx, group_state); - } - - // use NONE as default control - const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; - for (Phase phase : all) { - if (!group_state.has_injection_control(group.name(), phase)) { - group_state.injection_control(group.name(), phase, Group::InjectionCMode::NONE); - } - } - if (!group_state.has_production_control(group.name())) { - group_state.production_control(group.name(), Group::ProductionCMode::NONE); - } - - const auto& events = schedule[reportStepIdx].wellgroup_events(); - if (group.isInjectionGroup() - && events.hasEvent(group.name(), ScheduleEvents::GROUP_INJECTION_UPDATE)) { - - for (Phase phase : all) { - if (!group.hasInjectionControl(phase)) - continue; - - const auto& controls = group.injectionControls(phase, summaryState); - group_state.injection_control(group.name(), phase, controls.cmode); - } - } - - if (group.isProductionGroup() - && events.hasEvent(group.name(), ScheduleEvents::GROUP_PRODUCTION_UPDATE)) { - const auto controls = group.productionControls(summaryState); - group_state.production_control(group.name(), controls.cmode); - } - - if (group.has_gpmaint_control(Group::ProductionCMode::RESV)) { - group_state.production_control(group.name(), Group::ProductionCMode::RESV); - } - for (Phase phase : all) { - if (group.has_gpmaint_control(phase, Group::InjectionCMode::RATE)) { - group_state.injection_control(group.name(), phase, Group::InjectionCMode::RATE); - } else if (group.has_gpmaint_control(phase, Group::InjectionCMode::RESV)) { - group_state.injection_control(group.name(), phase, Group::InjectionCMode::RESV); - } - } - - if (schedule[reportStepIdx].gconsale().has(group.name())) { - group_state.injection_control(group.name(), Phase::GAS, Group::InjectionCMode::SALE); - } -} - -template -void WellGroupHelpers:: -accumulateGroupEfficiencyFactor(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - Scalar& factor) -{ - factor *= group.getGroupEfficiencyFactor(); - if (group.parent() != "FIELD" && !group.parent().empty()) - accumulateGroupEfficiencyFactor( - schedule.getGroup(group.parent(), reportStepIdx), schedule, reportStepIdx, factor); -} - -template -Scalar -WellGroupHelpers:: -sumWellSurfaceRates(const Group& group, - const Schedule& schedule, - const WellState& wellState, - const int reportStepIdx, - const int phasePos, - const bool injector, - const SummaryState& summaryState) -{ - return sumWellPhaseRates(false, group, schedule, wellState, summaryState, reportStepIdx, phasePos, injector); -} - -template -Scalar -WellGroupHelpers:: -sumWellResRates(const Group& group, - const Schedule& schedule, - const WellState& wellState, - const int reportStepIdx, - const int phasePos, - const bool injector, - const SummaryState& summaryState) -{ - return sumWellPhaseRates(true, group, schedule, wellState, summaryState, reportStepIdx, phasePos, injector); -} - -template -Scalar -WellGroupHelpers:: -sumSolventRates(const Group& group, - const Schedule& schedule, - const WellState& wellState, - const int reportStepIdx, - const bool injector) -{ - Scalar rate = 0.0; - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - const auto& gefac = groupTmp.getGroupEfficiencyFactor(); - rate += gefac * sumSolventRates(groupTmp, schedule, wellState, reportStepIdx, injector); - } - - for (const std::string& wellName : group.wells()) { - const auto well_index = wellState.index(wellName); - if (!well_index.has_value()) - continue; - - if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once - { - continue; - } - - const auto& wellEcl = schedule.getWell(wellName, reportStepIdx); - // only count producers or injectors - if ((wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector)) - continue; - - const auto& ws = wellState.well(well_index.value()); - - if (ws.status == Well::Status::SHUT) - continue; - - const Scalar factor = wellEcl.getEfficiencyFactor() * - wellState[wellEcl.name()].efficiency_scaling_factor; - if (injector) - rate += factor * ws.sum_solvent_rates(); - else - rate -= factor * ws.sum_solvent_rates(); - } - return rate; -} - -template -void WellGroupHelpers:: -updateGroupTargetReduction(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const bool isInjector, - const GuideRate& guide_rate, - const WellState& wellState, - const SummaryState& summaryState, - GroupState& group_state, - std::vector& groupTargetReduction) -{ - OPM_TIMEFUNCTION(); - const int np = wellState.numPhases(); - for (const std::string& subGroupName : group.groups()) { - std::vector subGroupTargetReduction(np, 0.0); - const Group& subGroup = schedule.getGroup(subGroupName, reportStepIdx); - updateGroupTargetReduction(subGroup, - schedule, - reportStepIdx, - isInjector, - guide_rate, - wellState, - summaryState, - group_state, - subGroupTargetReduction); - - const Scalar subGroupEfficiency = subGroup.getGroupEfficiencyFactor(); - - // accumulate group contribution from sub group - if (isInjector) { - const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; - for (Phase phase : all) { - - const auto& pu = wellState.phaseUsageInfo(); - int phase_pos = -1; - switch (phase) { - case Phase::GAS: - if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { - phase_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); - } - break; - case Phase::OIL: - if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { - phase_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); - } - break; - case Phase::WATER: - if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { - phase_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); - } - break; - default: - // just to avoid warning - throw std::invalid_argument("unhandled phase enum"); - } - - // the phase is not present - if (phase_pos == -1) - continue; - - const Group::InjectionCMode& currentGroupControl - = group_state.injection_control(subGroup.name(), phase); - const bool individual_control = (currentGroupControl != Group::InjectionCMode::FLD - && currentGroupControl != Group::InjectionCMode::NONE); - const int num_group_controlled_wells - = groupControlledWells(schedule, wellState, group_state, reportStepIdx, subGroupName, "", !isInjector, phase); - if (individual_control || num_group_controlled_wells == 0) { - groupTargetReduction[phase_pos] - += subGroupEfficiency * sumWellSurfaceRates(subGroup, schedule, wellState, reportStepIdx, phase_pos, isInjector, summaryState); - } else { - // Accumulate from this subgroup only if no group guide rate is set for it. - if (!guide_rate.has(subGroupName, phase)) { - groupTargetReduction[phase_pos] += subGroupEfficiency * subGroupTargetReduction[phase_pos]; - } - } - } - } else { - const Group::ProductionCMode& currentGroupControl = group_state.production_control(subGroupName); - const bool individual_control = (currentGroupControl != Group::ProductionCMode::FLD - && currentGroupControl != Group::ProductionCMode::NONE); - const int num_group_controlled_wells - = groupControlledWells(schedule, wellState, group_state, reportStepIdx, subGroupName, "", !isInjector, /*injectionPhaseNotUsed*/Phase::OIL); - if (individual_control || num_group_controlled_wells == 0) { - for (int phase = 0; phase < np; phase++) { - groupTargetReduction[phase] - += subGroupEfficiency * sumWellSurfaceRates(subGroup, schedule, wellState, reportStepIdx, phase, isInjector, summaryState); - } - } else { - // The subgroup may participate in group control. - if (!guide_rate.has(subGroupName)) { - // Accumulate from this subgroup only if no group guide rate is set for it. - for (int phase = 0; phase < np; phase++) { - groupTargetReduction[phase] += subGroupEfficiency * subGroupTargetReduction[phase]; - } - } - } - } - } - - for (const std::string& wellName : group.wells()) { - const auto& wellTmp = schedule.getWell(wellName, reportStepIdx); - - if (wellTmp.isProducer() && isInjector) - continue; - - if (wellTmp.isInjector() && !isInjector) - continue; - - const auto well_index = wellState.index(wellName); - if (!well_index.has_value()) - continue; - - if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once - { - continue; - } - - const auto& ws = wellState.well(well_index.value()); - if (ws.status == Well::Status::SHUT) - continue; - - const Scalar efficiency = wellTmp.getEfficiencyFactor() * - wellState[wellTmp.name()].efficiency_scaling_factor; - - // add contribution from wells not under group control - if (isInjector) { - if (ws.injection_cmode != Well::InjectorCMode::GRUP) - for (int phase = 0; phase < np; phase++) { - groupTargetReduction[phase] += ws.surface_rates[phase] * efficiency; - } - } else { - if ((ws.production_cmode != Well::ProducerCMode::GRUP)){ - if (!group.as_choke()) { - for (int phase = 0; phase < np; phase++) { - groupTargetReduction[phase] -= ws.surface_rates[phase] * efficiency; - } - } - } - } - } - if (isInjector) - group_state.update_injection_reduction_rates(group.name(), groupTargetReduction); - else - group_state.update_production_reduction_rates(group.name(), groupTargetReduction); -} - -template -void WellGroupHelpers:: -updateWellRatesFromGroupTargetScale(const Scalar scale, - const Group& group, - const Schedule& schedule, - const int reportStepIdx, - bool isInjector, - const GroupState& group_state, - WellState& wellState) -{ - OPM_TIMEFUNCTION(); - for (const std::string& groupName : group.groups()) { - bool individual_control = false; - if (isInjector) { - const Phase all[] = {Phase::WATER, Phase::OIL, Phase::GAS}; - for (Phase phase : all) { - const Group::InjectionCMode& currentGroupControl - = group_state.injection_control(groupName, phase); - individual_control = individual_control || (currentGroupControl != Group::InjectionCMode::FLD - && currentGroupControl != Group::InjectionCMode::NONE); - } - } else { - const Group::ProductionCMode& currentGroupControl = group_state.production_control(groupName); - individual_control = (currentGroupControl != Group::ProductionCMode::FLD - && currentGroupControl != Group::ProductionCMode::NONE); - } - if (!individual_control) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateWellRatesFromGroupTargetScale(scale, groupTmp, schedule, reportStepIdx, isInjector, group_state, wellState); - } - } - - const int np = wellState.numPhases(); - for (const std::string& wellName : group.wells()) { - const auto& wellTmp = schedule.getWell(wellName, reportStepIdx); - - if (wellTmp.isProducer() && isInjector) - continue; - - if (wellTmp.isInjector() && !isInjector) - continue; - - const auto well_index = wellState.index(wellName); - if (!well_index.has_value()) - continue; - - auto& ws = wellState.well(well_index.value()); - if (ws.status == Well::Status::SHUT) - continue; - - // scale rates - if (isInjector) { - if (ws.injection_cmode == Well::InjectorCMode::GRUP) - for (int phase = 0; phase < np; phase++) { - ws.surface_rates[phase] *= scale; - } - } else { - if (ws.production_cmode == Well::ProducerCMode::GRUP) { - for (int phase = 0; phase < np; phase++) { - ws.surface_rates[phase] *= scale; - } - } - } - } - - -} - -template -void WellGroupHelpers:: -updateVREPForGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellState& wellState, - GroupState& group_state, - const SummaryState& summaryState) -{ - OPM_TIMEFUNCTION(); - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateVREPForGroups(groupTmp, schedule, reportStepIdx, wellState, group_state, summaryState); - } - const int np = wellState.numPhases(); - Scalar resv = 0.0; - for (int phase = 0; phase < np; ++phase) { - resv += sumWellPhaseRates(true, - group, - schedule, - wellState, - summaryState, - reportStepIdx, - phase, - /*isInjector*/ false); - } - group_state.update_injection_vrep_rate(group.name(), resv); -} - -template -void WellGroupHelpers:: -updateReservoirRatesInjectionGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellState& wellState, - GroupState& group_state, - const SummaryState& summaryState) -{ - OPM_TIMEFUNCTION(); - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateReservoirRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellState, group_state, summaryState); - } - const int np = wellState.numPhases(); - std::vector resv(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - resv[phase] = sumWellPhaseRates(true, - group, - schedule, - wellState, - summaryState, - reportStepIdx, - phase, - /*isInjector*/ true); - } - group_state.update_injection_reservoir_rates(group.name(), resv); -} - -template -void WellGroupHelpers:: -updateSurfaceRatesInjectionGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellState& wellState, - GroupState& group_state, - const SummaryState& summaryState) -{ - OPM_TIMEFUNCTION(); - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateSurfaceRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellState, group_state, summaryState); - } - const int np = wellState.numPhases(); - std::vector rates(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - rates[phase] = sumWellPhaseRates(false, - group, - schedule, - wellState, - summaryState, - reportStepIdx, - phase, - /*isInjector*/ true); - } - group_state.update_injection_surface_rates(group.name(), rates); -} - -template -void WellGroupHelpers:: -updateWellRates(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellState& wellStateNupcol, - WellState& wellState) -{ - OPM_TIMEFUNCTION(); - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateWellRates(groupTmp, schedule, reportStepIdx, wellStateNupcol, wellState); - } - const int np = wellState.numPhases(); - for (const std::string& wellName : group.wells()) { - std::vector rates(np, 0.0); - const auto well_index = wellState.index(wellName); - if (well_index.has_value()) { // the well is found on this node - const auto& wellTmp = schedule.getWell(wellName, reportStepIdx); - int sign = 1; - // production wellRates are negative. The users of currentWellRates uses the convention in - // opm-common that production and injection rates are positive. - if (!wellTmp.isInjector()) - sign = -1; - const auto& ws = wellStateNupcol.well(well_index.value()); - for (int phase = 0; phase < np; ++phase) { - rates[phase] = sign * ws.surface_rates[phase]; - } - } - wellState.setCurrentWellRates(wellName, rates); - } -} - -template -void WellGroupHelpers:: -updateGroupProductionRates(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellState& wellState, - GroupState& group_state, - const SummaryState& summaryState) -{ - OPM_TIMEFUNCTION(); - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateGroupProductionRates(groupTmp, schedule, reportStepIdx, wellState, group_state, summaryState); - } - const int np = wellState.numPhases(); - std::vector rates(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, summaryState, reportStepIdx, phase, /*isInjector*/ false); - } - group_state.update_production_rates(group.name(), rates); -} - -template -void WellGroupHelpers:: -updateNetworkLeafNodeProductionRates(const Schedule& schedule, - const int reportStepIdx, - const WellState& wellState, - GroupState& group_state, - const SummaryState& summaryState) -{ - const auto& network = schedule[reportStepIdx].network(); - if (network.active()) { - const int np = wellState.numPhases(); - for (const auto& group_name : network.leaf_nodes()) { - std::vector network_rates(np, 0.0); - if (schedule[reportStepIdx].groups.has(group_name)) { // Allow empty leaf nodes that are not groups - const auto& group = schedule[reportStepIdx].groups.get(group_name); - if (group.numWells() > 0) { - for (int phase = 0; phase < np; ++phase) { - network_rates[phase] = sumWellPhaseRates(false, group, schedule, wellState, summaryState, reportStepIdx, phase, /*isInjector*/ false, /*network*/ true); - } - } - } - group_state.update_network_leaf_node_production_rates(group_name, network_rates); - } - } -} - -template -void WellGroupHelpers:: -updateREINForGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const SummaryState& st, - const WellState& wellState, - GroupState& group_state, - bool sum_rank) -{ - OPM_TIMEFUNCTION(); - const int np = wellState.numPhases(); - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateREINForGroups(groupTmp, schedule, reportStepIdx, st, wellState, group_state, sum_rank); - } - - std::vector rein(np, 0.0); - for (int phase = 0; phase < np; ++phase) { - rein[phase] = sumWellPhaseRates(false, group, schedule, wellState, st, reportStepIdx, phase, /*isInjector*/ false); - } - - // add import rate and subtract consumption rate for group for gas - if (sum_rank) { - const auto& pu = wellState.phaseUsageInfo(); - if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { - const auto& [consumption_rate, import_rate] = group_state.gconsump_rates(group.name()); - const auto ig = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); - rein[ig] += import_rate; - rein[ig] -= consumption_rate; - } - } - - group_state.update_injection_rein_rates(group.name(), rein); -} - -template -template -void WellGroupHelpers:: -updateGpMaintTargetForGroups(const Group& group, - const Schedule& schedule, - const RegionalValues& regional_values, - const int reportStepIdx, - const double dt, - const WellState& well_state, - GroupState& group_state) -{ - OPM_TIMEFUNCTION(); - for (const std::string& groupName : group.groups()) { - const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); - updateGpMaintTargetForGroups(groupTmp, schedule, regional_values, reportStepIdx, dt, well_state, group_state); - } - const auto& gpm = group.gpmaint(); - if (!gpm) - return; - - const auto& region = gpm->region(); - if (!region) - return; - - const auto [name, number] = *region; - const Scalar error = gpm->pressure_target() - regional_values.at(name)->pressure(number); - Scalar current_rate = 0.0; - const auto& pu = well_state.phaseUsageInfo(); - bool injection = true; - Scalar sign = 1.0; - switch (gpm->flow_target()) { - case GPMaint::FlowTarget::RESV_PROD: - { - current_rate = -group_state.injection_vrep_rate(group.name()); - injection = false; - sign = -1.0; - break; - } - case GPMaint::FlowTarget::RESV_OINJ: - { - if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { - const auto io = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); - current_rate = group_state.injection_reservoir_rates(group.name())[io]; - } - break; - } - case GPMaint::FlowTarget::RESV_WINJ: - { - if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { - const auto iw = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); - current_rate = group_state.injection_reservoir_rates(group.name())[iw]; - } - break; - } - case GPMaint::FlowTarget::RESV_GINJ: - { - if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { - const auto ig = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); - current_rate = group_state.injection_reservoir_rates(group.name())[ig]; - } - break; - } - case GPMaint::FlowTarget::SURF_OINJ: - { - if (pu.phaseIsActive(IndexTraits::oilPhaseIdx)) { - const auto io = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); - current_rate = group_state.injection_surface_rates(group.name())[io]; - } - break; - } - case GPMaint::FlowTarget::SURF_WINJ: - { - if (pu.phaseIsActive(IndexTraits::waterPhaseIdx)) { - const auto iw = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); - current_rate = group_state.injection_surface_rates(group.name())[iw]; - } - break; - } - case GPMaint::FlowTarget::SURF_GINJ: - { - if (pu.phaseIsActive(IndexTraits::gasPhaseIdx)) { - const auto ig = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); - current_rate = group_state.injection_surface_rates(group.name())[ig]; - } - break; - } - default: - throw std::invalid_argument("Invalid Flow target type in GPMAINT"); - } - auto& gpmaint_state = group_state.gpmaint(group.name()); - // we only activate gpmaint if pressure is lower than the target regional pressure for injectors - // (i.e. error > 0) and higher for producers. - bool activate = (injection && error > 0) || (!injection && error < 0); - Scalar rate = 0.0; - if (activate) { - rate = gpm->rate(gpmaint_state, current_rate, error, dt); - } else { - gpm->resetState(gpmaint_state); - } - group_state.update_gpmaint_target(group.name(), std::max(Scalar{0.0}, sign * rate)); -} - -template -std::map -WellGroupHelpers:: -computeNetworkPressures(const Network::ExtNetwork& network, - const WellStateType& well_state, - const GroupState& group_state, - const VFPProdProperties& vfp_prod_props, - const Schedule& schedule, - const Parallel::Communication& comm, - const int report_time_step) -{ - // 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 (!group_state.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 = schedule.getGroup(node, report_time_step); - Scalar alq = 0.0; - for (const std::string& wellname : group.wells()) { - const Well& well = schedule.getWell(wellname, report_time_step); - if (well.isInjector() || !well_state.isOpen(wellname)) continue; - const Scalar efficiency = well.getEfficiencyFactor(/*network*/ true) * well_state.getGlobalEfficiencyScalingFactor(wellname); - const auto& well_index = well_state.index(wellname); - if (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] = group_state.well_group_thp(node); - } else { - node_pressures[node] = up_press; - } - } - } - } - } - return node_pressures; -} - -template -GuideRate::RateVector -WellGroupHelpers:: -getWellRateVector(const WellState& well_state, - const std::string& name) -{ - return getGuideRateVector(well_state.currentWellRates(name), well_state.phaseUsageInfo()); -} - -template -GuideRate::RateVector -WellGroupHelpers:: -getProductionGroupRateVector(const GroupState& group_state, - const PhaseUsageInfo& pu, - const std::string& group_name) -{ - return getGuideRateVector(group_state.production_rates(group_name), pu); -} - -template -Scalar -WellGroupHelpers:: -getGuideRate(const std::string& name, - const Schedule& schedule, - const WellState& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const GuideRateModel::Target target) -{ - if (schedule.hasWell(name, reportStepIdx)) { - if (guideRate->has(name) || guideRate->hasPotentials(name)) { - return guideRate->get(name, target, getWellRateVector(wellState, name)); - } else { - return 0.0; - } - } - - if (guideRate->has(name)) { - return guideRate->get(name, target, getProductionGroupRateVector(group_state, wellState.phaseUsageInfo(), name)); - } - - Scalar totalGuideRate = 0.0; - const Group& group = schedule.getGroup(name, reportStepIdx); - - for (const std::string& groupName : group.groups()) { - const Group::ProductionCMode& currentGroupControl = group_state.production_control(groupName); - if (currentGroupControl == Group::ProductionCMode::FLD - || currentGroupControl == Group::ProductionCMode::NONE) { - // accumulate from sub wells/groups - totalGuideRate += getGuideRate(groupName, schedule, wellState, group_state, reportStepIdx, guideRate, target); - } - } - - for (const std::string& wellName : group.wells()) { - const auto& wellTmp = schedule.getWell(wellName, reportStepIdx); - - if (wellTmp.isInjector()) - continue; - - const auto well_index = wellState.index(wellName); - if (!well_index.has_value()) - continue; - - if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once - { - continue; - } - - const auto& ws = wellState.well(well_index.value()); - if (ws.status == Well::Status::SHUT) - continue; - - // Only count wells under group control or the ru - if (!wellState.isProductionGrup(wellName)) - continue; - - totalGuideRate += getGuideRate(wellName, schedule, wellState, group_state, - reportStepIdx, guideRate, target); - - } - return totalGuideRate; -} - -template -Scalar -WellGroupHelpers:: -getGuideRateInj(const std::string& name, - const Schedule& schedule, - const WellState& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const GuideRateModel::Target target, - const Phase& injectionPhase) -{ - if (schedule.hasWell(name, reportStepIdx)) { - return getGuideRate(name, schedule, wellState, group_state, - reportStepIdx, guideRate, target); - } - - if (guideRate->has(name, injectionPhase)) { - return guideRate->get(name, injectionPhase); - } - - Scalar totalGuideRate = 0.0; - const Group& group = schedule.getGroup(name, reportStepIdx); - - for (const std::string& groupName : group.groups()) { - const Group::InjectionCMode& currentGroupControl - = group_state.injection_control(groupName, injectionPhase); - if (currentGroupControl == Group::InjectionCMode::FLD - || currentGroupControl == Group::InjectionCMode::NONE) { - // accumulate from sub wells/groups - totalGuideRate += getGuideRateInj(groupName, schedule, wellState, group_state, reportStepIdx, guideRate, target, injectionPhase); - } - } - - for (const std::string& wellName : group.wells()) { - const auto& wellTmp = schedule.getWell(wellName, reportStepIdx); - - if (!wellTmp.isInjector()) - continue; - - const auto well_index = wellState.index(wellName); - if (!well_index.has_value()) - continue; - - if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once - { - continue; - } - - const auto& ws = wellState.well(well_index.value()); - if (ws.status == Well::Status::SHUT) - continue; - - // Only count wells under group control or the ru - if (!wellState.isInjectionGrup(wellName)) - continue; - - totalGuideRate += guideRate->get(wellName, target, getWellRateVector(wellState, wellName)); - } - return totalGuideRate; -} - -template -int WellGroupHelpers:: -updateGroupControlledWells(const Schedule& schedule, - const WellStateType& well_state, - GroupState& group_state, - const SummaryState& summary_state, - const GuideRate* guideRate, - const int report_step, - const std::string& group_name, - const bool is_production_group, - const Phase injection_phase) -{ - OPM_TIMEFUNCTION(); - const Group& group = schedule.getGroup(group_name, report_step); - int num_wells = 0; - for (const std::string& child_group : group.groups()) { - - bool included = false; - if (is_production_group) { - const auto ctrl = group_state.production_control(child_group); - included = included || (ctrl == Group::ProductionCMode::FLD) || (ctrl == Group::ProductionCMode::NONE); - } else { - const auto ctrl = group_state.injection_control(child_group, injection_phase); - included = included || (ctrl == Group::InjectionCMode::FLD) || (ctrl == Group::InjectionCMode::NONE); - } - - if (included) { - num_wells - += updateGroupControlledWells(schedule, well_state, group_state, summary_state, guideRate, report_step, child_group, is_production_group, injection_phase); - } else { - updateGroupControlledWells(schedule, well_state, group_state, summary_state, guideRate, report_step, child_group, is_production_group, injection_phase); - } - } - for (const std::string& child_well : group.wells()) { - bool included = false; - const Well& well = schedule.getWell(child_well, report_step); - if (is_production_group && well.isProducer()) { - included = included || well_state.isProductionGrup(child_well) || group.as_choke(); - } else if (!is_production_group && !well.isProducer()) { - const auto& well_controls = well.injectionControls(summary_state); - auto injectorType = well_controls.injector_type; - if ( (injection_phase == Phase::WATER && injectorType == InjectorType::WATER ) || - (injection_phase == Phase::OIL && injectorType == InjectorType::OIL ) || - (injection_phase == Phase::GAS && injectorType == InjectorType::GAS )) - { - included = included || well_state.isInjectionGrup(child_well); - } - } - const auto ctrl1 = group_state.production_control(group.name()); - if (group.as_choke() && ((ctrl1 == Group::ProductionCMode::FLD) || (ctrl1 == Group::ProductionCMode::NONE))){ - // The auto choke group has not own group control but inherits control from an ancestor group. - // Number of wells should be calculated as zero when wells of auto choke group do not deliver target. - // This behaviour is then similar to no-autochoke group with wells not on GRUP control. - // The rates of these wells are summed up. The parent group target is reduced with this rate. - // This reduced target becomes the target of the other child group of this parent. - const auto num_phases = well_state.numPhases(); - std::vector rates(num_phases, 0.0); - for (int phase_pos = 0; phase_pos < num_phases; ++phase_pos) { - rates[phase_pos] = WellGroupHelpers::sumWellSurfaceRates(group, - schedule, - well_state, - report_step, - phase_pos, - false, - summary_state); - } - - // Get the ancestor of the auto choke group that has group control (cmode != FLD, NONE) - const auto& control_group_name = control_group(group, group_state, report_step, schedule); - const auto& control_group = schedule.getGroup(control_group_name, report_step); - const auto& ctrl = control_group.productionControls(summary_state); - const auto& control_group_cmode = ctrl.cmode; - - const auto& group_guide_rate = group.productionControls(summary_state).guide_rate; - - if (group_guide_rate > 0) { - // Guide rate is not default for the auto choke group - Scalar gratTargetFromSales = 0.0; - if (group_state.has_grat_sales_target(control_group_name)) - gratTargetFromSales = group_state.grat_sales_target(control_group_name); - - std::vector resv_coeff(num_phases, 1.0); - WGHelpers::TargetCalculator tcalc(control_group_cmode, - well_state.phaseUsageInfo(), - resv_coeff, - gratTargetFromSales, - group.name(), - group_state, - group.has_gpmaint_control(control_group_cmode)); - auto deferred_logger = Opm::DeferredLogger(); - const auto& control_group_target = tcalc.groupTarget(ctrl, deferred_logger); - - // Calculates the guide rate of the parent group with control. - // It is allowed that the guide rate of this group is defaulted. The guide rate will be derived from the children groups - const auto& control_group_guide_rate = getGuideRate(control_group_name, - schedule, - well_state, - group_state, - report_step, - guideRate, - tcalc.guideTargetMode()); - - if (control_group_guide_rate > 0) { - // Target rate for the auto choke group - const Scalar target_rate = control_group_target * group_guide_rate / control_group_guide_rate; - const Scalar current_rate = tcalc.calcModeRateFromRates(rates); - - if (current_rate < target_rate) - included = false; - } - } - } - - if (included) { - ++num_wells; - } - } - if (is_production_group) { - group_state.update_number_of_wells_under_group_control(group_name, num_wells); - } else { - group_state.update_number_of_wells_under_inj_group_control(group_name, injection_phase, num_wells); - } - - return num_wells; -} - - -template -int WellGroupHelpers:: -groupControlledWells(const Schedule& schedule, - const WellStateType& well_state, - const GroupState& group_state, - const int report_step, - const std::string& group_name, - const std::string& always_included_child, - const bool is_production_group, - const Phase injection_phase) -{ - auto num_wells = is_production_group ? group_state.number_of_wells_under_group_control(group_name) - : group_state.number_of_wells_under_inj_group_control(group_name, injection_phase); - if (schedule.hasWell(always_included_child, report_step)) { - const bool isInGroup = isInGroupChainTopBot(always_included_child, group_name, schedule, report_step); - const bool already_included = is_production_group ? well_state.isProductionGrup(always_included_child) - : well_state.isInjectionGrup(always_included_child); - if (!already_included && isInGroup) { - num_wells++; - } - } - return num_wells; -} - -template -std::vector -WellGroupHelpers:: -groupChainTopBot(const std::string& bottom, - const std::string& top, - const Schedule& schedule, - const int report_step) -{ - // Get initial parent, 'bottom' can be a well or a group. - std::string parent; - if (schedule.hasWell(bottom, report_step)) { - parent = schedule.getWell(bottom, report_step).groupName(); - } else { - parent = schedule.getGroup(bottom, report_step).parent(); - } - - // Build the chain from bottom to top. - std::vector chain; - chain.push_back(bottom); - chain.push_back(parent); - while (parent != top) { - parent = schedule.getGroup(parent, report_step).parent(); - chain.push_back(parent); - } - assert(chain.back() == top); - - // Reverse order and return. - std::reverse(chain.begin(), chain.end()); - return chain; -} - -template -bool -WellGroupHelpers:: -isInGroupChainTopBot(const std::string& bottom, - const std::string& top, - const Schedule& schedule, - const int report_step) -{ - // Get initial parent, 'bottom' can be a well or a group. - std::string parent; - if (schedule.hasWell(bottom, report_step)) { - parent = schedule.getWell(bottom, report_step).groupName(); - } else { - parent = schedule.getGroup(bottom, report_step).parent(); - } - - while (parent != top) { - parent = schedule.getGroup(parent, report_step).parent(); - if (parent == top) { - return true; - } else if (parent == "FIELD") { - return false; - } - } - return true; -} - -template -std::string -WellGroupHelpers:: -control_group(const Group& group, - const GroupState& group_state, - const int reportStepIdx, - const Schedule& schedule) -{ - const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name()); - - if (currentGroupControl == Group::ProductionCMode::FLD || currentGroupControl == Group::ProductionCMode::NONE) { - const auto& parent_name = group.control_group(); - if (parent_name) { - const auto& parent = schedule.getGroup(parent_name.value(), reportStepIdx); - return control_group(parent, - group_state, - reportStepIdx, - schedule); - } - } - - return group.name(); -} - -template -std::pair -WellGroupHelpers:: -checkGroupConstraintsProd(const std::string& name, - const std::string& parent, - const Group& group, - const WellState& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - const bool check_guide_rate, - DeferredLogger& deferred_logger) -{ - // When called for a well ('name' is a well name), 'parent' - // will be the name of 'group'. But if we recurse, 'name' and - // 'parent' will stay fixed while 'group' will be higher up - // in the group tree. - // efficiencyfactor is the well efficiency factor for the first group the well is - // part of. Later it is the accumulated factor including the group efficiency factor - // of the child of group. - - // NOTE: if name is a group, the rates-argument is the portion of the 'full' group rates that - // is potentially available for control by an ancestor, i.e. full rates minus reduction rates - OPM_TIMEFUNCTION(); - const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name()); - - if (currentGroupControl == Group::ProductionCMode::FLD || currentGroupControl == Group::ProductionCMode::NONE) { - // Return if we are not available for parent group. - if (!group.productionGroupControlAvailable()) { - return std::make_pair(false, 1); - } - // Otherwise: check production share of parent's control. - const auto& parentGroup = schedule.getGroup(group.parent(), reportStepIdx); - return checkGroupConstraintsProd(name, - parent, - parentGroup, - wellState, - group_state, - reportStepIdx, - guideRate, - rates, - efficiencyFactor * group.getGroupEfficiencyFactor(), - schedule, - summaryState, - resv_coeff, - check_guide_rate, - deferred_logger); - } - - // This can be false for FLD-controlled groups, we must therefore - // check for FLD first (done above). - if (!group.isProductionGroup()) { - return std::make_pair(false, Scalar(1.0)); - } - - // If we are here, we are at the topmost group to be visited in the recursion. - // This is the group containing the control we will check against. - - // gconsale may adjust the grat target. - // the adjusted rates is send to the targetCalculator - Scalar gratTargetFromSales = 0.0; - if (group_state.has_grat_sales_target(group.name())) - gratTargetFromSales = group_state.grat_sales_target(group.name()); - - WGHelpers::TargetCalculator tcalc(currentGroupControl, - wellState.phaseUsageInfo(), - resv_coeff, - gratTargetFromSales, - group.name(), - group_state, - group.has_gpmaint_control(currentGroupControl)); - - WGHelpers::FractionCalculator fcalc(schedule, - wellState, - group_state, - summaryState, - reportStepIdx, - guideRate, - tcalc.guideTargetMode(), - true, - Phase::OIL); - - auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); }; - - auto localReduction = [&](const std::string& group_name) - { - const std::vector& groupTargetReductions = - group_state.production_reduction_rates(group_name); - return tcalc.calcModeRateFromRates(groupTargetReductions); - }; - - auto localCurrentRate = [&](const std::string& group_name) - { - const std::vector& groupSurfaceRates = - group_state.production_rates(group_name); - return tcalc.calcModeRateFromRates(groupSurfaceRates); - }; - - std::optional ctrl; - if (!group.has_gpmaint_control(currentGroupControl)) - ctrl = group.productionControls(summaryState); - - const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); - // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. - // Then ... - // TODO finish explanation. - const Scalar current_rate_available - = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. - const auto chain = groupChainTopBot(name, group.name(), schedule, reportStepIdx); - // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. - const std::size_t num_ancestors = chain.size() - 1; - - // check whether guide rate is violated - if (check_guide_rate) { - for (std::size_t ii = 1; ii < num_ancestors; ++ii) { - if (guideRate->has(chain[ii])) { - const auto& guided_group = chain[ii]; - const Scalar grefficiency - = schedule.getGroup(guided_group, reportStepIdx).getGroupEfficiencyFactor(); - const Scalar currentRateFraction = grefficiency * localCurrentRate(guided_group) / (localCurrentRate(chain[ii-1])); - const Scalar guiderateFraction = localFraction(guided_group); - // we add a factor here to avoid switching due to numerical instability - const Scalar factor = 1.01; - if (currentRateFraction > (guiderateFraction * factor)) { - return std::make_pair(true, guiderateFraction / currentRateFraction); - } - } - } - } - - if (schedule.hasWell(name) && wellState.well(name).group_target) { // for wells we already have computed the target - const Scalar current_well_rate_available = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. - const Scalar group_target_rate_available = *wellState.well(name).group_target; - Scalar scale = 1.0; - if (current_well_rate_available > 1e-12) { - scale = group_target_rate_available / current_well_rate_available; - } - - return std::make_pair(current_well_rate_available > group_target_rate_available, scale); - } - - // we need to find out the level where the current well is applied to the local reduction - std::size_t local_reduction_level = 0; - for (std::size_t ii = 1; ii < num_ancestors; ++ii) { - const int num_gr_ctrl = groupControlledWells(schedule, - wellState, - group_state, - reportStepIdx, - chain[ii], - "", - /*is_producer*/ true, - /*injectionPhaseNotUsed*/ Phase::OIL); - if (guideRate->has(chain[ii]) && num_gr_ctrl > 0) { - local_reduction_level = ii; - } - } - - // Compute portion of target corresponding to current_rate_available - Scalar target = orig_target; - for (std::size_t ii = 0; ii < num_ancestors; ++ii) { - if ((ii == 0) || guideRate->has(chain[ii])) { - // Apply local reductions only at the control level - // (top) and for levels where we have a specified - // group guide rate. - if (local_reduction_level >= ii ) { - target -= localReduction(chain[ii]); - } - // Add my reduction back at the level where it is included in the local reduction - if (local_reduction_level == ii ) { - target += current_rate_available * efficiencyFactor; - } - } - target *= localFraction(chain[ii + 1]); - } - // Avoid negative target rates comming from too large local reductions. - const Scalar target_rate_available = std::max(Scalar(1e-12), target / efficiencyFactor); - - Scalar scale = 1.0; - if (current_rate_available > 1e-12) - scale = target_rate_available / current_rate_available; - - return std::make_pair(current_rate_available > target_rate_available, scale); -} - - -template -Scalar -WellGroupHelpers:: -getWellGroupTargetProducer(const std::string& name, - const std::string& parent, - const Group& group, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - DeferredLogger& deferred_logger) -{ - // This function computes a wells group target. - // 'parent' will be the name of 'group'. But if we recurse, 'name' and - // 'parent' will stay fixed while 'group' will be higher up - // in the group tree. - // Eficiencyfactor is the well efficiency factor for the first group the well is - // part of. Later it is the accumulated factor including the group efficiency factor - // of the child of group. - OPM_TIMEFUNCTION(); - const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name()); - - if (currentGroupControl == Group::ProductionCMode::FLD || currentGroupControl == Group::ProductionCMode::NONE) { - // Return if we are not available for parent group. - if (!group.productionGroupControlAvailable()) { - return std::numeric_limits::max(); - } - // Otherwise: check production share of parent's control. - const auto& parentGroup = schedule.getGroup(group.parent(), reportStepIdx); - return getWellGroupTargetProducer(name, - parent, - parentGroup, - wellState, - group_state, - reportStepIdx, - guideRate, - rates, - efficiencyFactor * group.getGroupEfficiencyFactor(), - schedule, - summaryState, - resv_coeff, - deferred_logger); - } - - // This can be false for FLD-controlled groups, we must therefore - // check for FLD first (done above). - if (!group.isProductionGroup()) { - return std::numeric_limits::max(); - } - - // If we are here, we are at the topmost group to be visited in the recursion. - // This is the group containing the control we will check against. - - // gconsale may adjust the grat target. - // the adjusted rates is send to the targetCalculator - Scalar gratTargetFromSales = 0.0; - if (group_state.has_grat_sales_target(group.name())) - gratTargetFromSales = group_state.grat_sales_target(group.name()); - - WGHelpers::TargetCalculator tcalc(currentGroupControl, - wellState.phaseUsageInfo(), - resv_coeff, - gratTargetFromSales, - group.name(), - group_state, - group.has_gpmaint_control(currentGroupControl)); - - WGHelpers::FractionCalculator fcalc(schedule, - wellState, - group_state, - summaryState, - reportStepIdx, - guideRate, - tcalc.guideTargetMode(), - true, - Phase::OIL); - auto localFraction = [&](const std::string& child, const std::string& always_incluced_name) { return fcalc.localFraction(child, always_incluced_name); }; - - auto localReduction = [&](const std::string& group_name) - { - const std::vector& groupTargetReductions = - group_state.production_reduction_rates(group_name); - return tcalc.calcModeRateFromRates(groupTargetReductions); - }; - - std::optional ctrl; - if (!group.has_gpmaint_control(currentGroupControl)) - ctrl = group.productionControls(summaryState); - - const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); - // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. - // Then ... - // TODO finish explanation. - const Scalar current_rate_available - = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. - const auto chain = groupChainTopBot(name, group.name(), schedule, reportStepIdx); - // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. - const std::size_t num_ancestors = chain.size() - 1; - // we need to find out the level where the local reduction is applied - std::size_t local_reduction_level = 0; - for (std::size_t ii = 1; ii < num_ancestors; ++ii) { - const int num_gr_ctrl = groupControlledWells(schedule, - wellState, - group_state, - reportStepIdx, - chain[ii], - "", - /*is_producer*/ true, - /*injectionPhaseNotUsed*/ Phase::OIL); - if (guideRate->has(chain[ii]) && num_gr_ctrl > 0) { - local_reduction_level = ii; - } - } - // Compute portion of target corresponding to current_rate_available - Scalar target = orig_target; - for (std::size_t ii = 0; ii < num_ancestors; ++ii) { - if ((ii == 0) || guideRate->has(chain[ii])) { - // Apply local reductions only at the control level - // (top) and for levels where we have a specified - // group guide rate. - if (local_reduction_level >= ii) { - target -= localReduction(chain[ii]); - } - // If we are under individual control we need to add the wells rate back at the level where it is included in the local reduction - if (local_reduction_level == ii && !wellState.isProductionGrup(name)) { - target += current_rate_available * efficiencyFactor; - } - } - if (wellState.isProductionGrup(name)) { - target *= localFraction(chain[ii + 1], chain[ii + 1]); - } else { - target *= localFraction(chain[ii + 1], name); - } - } - // Avoid negative target rates comming from too large local reductions. - return std::max(Scalar(0.0), target / efficiencyFactor);; -} - -template -std::pair -WellGroupHelpers:: -checkGroupConstraintsInj(const std::string& name, - const std::string& parent, - const Group& group, - const WellState& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - Phase injectionPhase, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - const bool check_guide_rate, - DeferredLogger& deferred_logger) -{ - // When called for a well ('name' is a well name), 'parent' - // will be the name of 'group'. But if we recurse, 'name' and - // 'parent' will stay fixed while 'group' will be higher up - // in the group tree. - // efficiencyfactor is the well efficiency factor for the first group the well is - // part of. Later it is the accumulated factor including the group efficiency factor - // of the child of group. - - // NOTE: if name is a group, the rates-argument is the portion of the 'full' group rates that - // is potentially available for control by an ancestor, i.e. full rates minus reduction rates - - auto currentGroupControl = group_state.injection_control(group.name(), injectionPhase); - - if (currentGroupControl == Group::InjectionCMode::FLD || currentGroupControl == Group::InjectionCMode::NONE) { - // Return if we are not available for parent group. - if (!group.injectionGroupControlAvailable(injectionPhase)) { - return std::make_pair(false, Scalar(1)); - } - // Otherwise: check production share of parent's control. - const auto& parentGroup = schedule.getGroup(group.parent(), reportStepIdx); - return checkGroupConstraintsInj(name, - parent, - parentGroup, - wellState, - group_state, - reportStepIdx, - guideRate, - rates, - injectionPhase, - efficiencyFactor * group.getGroupEfficiencyFactor(), - schedule, - summaryState, - resv_coeff, - check_guide_rate, - deferred_logger); - } - - // This can be false for FLD-controlled groups, we must therefore - // check for FLD first (done above). - if (!group.isInjectionGroup()) { - return std::make_pair(false, Scalar(1.0)); - } - - // If we are here, we are at the topmost group to be visited in the recursion. - // This is the group containing the control we will check against. - Scalar sales_target = 0; - if (schedule[reportStepIdx].gconsale().has(group.name())) { - const auto& gconsale = schedule[reportStepIdx].gconsale().get(group.name(), summaryState); - sales_target = gconsale.sales_target; - } - WGHelpers::InjectionTargetCalculator tcalc(currentGroupControl, - wellState.phaseUsageInfo(), - resv_coeff, - group.name(), - sales_target, - group_state, - injectionPhase, - group.has_gpmaint_control(injectionPhase, - currentGroupControl), - deferred_logger); - - WGHelpers::FractionCalculator fcalc(schedule, - wellState, - group_state, - summaryState, - reportStepIdx, - guideRate, - tcalc.guideTargetMode(), - false, - injectionPhase); - - auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); }; - - auto localReduction = [&](const std::string& group_name) - { - const std::vector& groupTargetReductions = - group_state.injection_reduction_rates(group_name); - return tcalc.calcModeRateFromRates(groupTargetReductions); - }; - - auto localCurrentRate = [&](const std::string& group_name) - { - const std::vector& groupSurfaceRates = - group_state.injection_surface_rates(group_name); - return tcalc.calcModeRateFromRates(groupSurfaceRates); - }; - - const auto chain = groupChainTopBot(name, group.name(), schedule, reportStepIdx); - // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. - const std::size_t num_ancestors = chain.size() - 1; - // we need to find out the level where the current well is applied to the local reduction - std::size_t local_reduction_level = 0; - for (std::size_t ii = 1; ii < num_ancestors; ++ii) { - const int num_gr_ctrl = groupControlledWells(schedule, - wellState, - group_state, - reportStepIdx, - chain[ii], - "", - /*is_producer*/ false, - injectionPhase); - if (guideRate->has(chain[ii], injectionPhase) && num_gr_ctrl > 0) { - local_reduction_level = ii; - } - } - - // check whether guide rate is violated - if (check_guide_rate) { - for (std::size_t ii = 1; ii < num_ancestors; ++ii) { - if (guideRate->has(chain[ii], injectionPhase)) { - const auto& guided_group = chain[ii]; - const Scalar grefficiency - = schedule.getGroup(guided_group, reportStepIdx).getGroupEfficiencyFactor(); - const Scalar currentRateFraction = grefficiency * localCurrentRate(guided_group) / - localCurrentRate(chain[ii-1]); - const Scalar guiderateFraction = localFraction(guided_group); - // we add a factor here to avoid switching due to numerical instability - const Scalar factor = 1.01; - if (currentRateFraction > (guiderateFraction * factor)) { - return std::make_pair(true, guiderateFraction / currentRateFraction); - } - } - } - } - - if (schedule.hasWell(name) && wellState.well(name).group_target) { // for wells we already have computed the target - Scalar scale = 1.0; - const Scalar group_target_rate_available = *wellState.well(name).group_target; - const Scalar current_well_rate_available = tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. - if (current_well_rate_available > 1e-12) { - scale = group_target_rate_available / current_well_rate_available; - } - return std::make_pair(current_well_rate_available > group_target_rate_available, scale); - } - - std::optional ctrl; - if (!group.has_gpmaint_control(injectionPhase, currentGroupControl)) - ctrl = group.injectionControls(injectionPhase, summaryState); - - - const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); - // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. - // Then ... - // TODO finish explanation. - const Scalar current_rate_available - = tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. - - // Compute portion of target corresponding to current_rate_available - Scalar target = orig_target; - for (std::size_t ii = 0; ii < num_ancestors; ++ii) { - if ((ii == 0) || guideRate->has(chain[ii], injectionPhase)) { - // Apply local reductions only at the control level - // (top) and for levels where we have a specified - // group guide rate. - if (local_reduction_level >= ii ) { - target -= localReduction(chain[ii]); - } - - // Add my reduction back at the level where it is included in the local reduction - if (local_reduction_level == ii ) { - target += current_rate_available * efficiencyFactor; - } - } - target *= localFraction(chain[ii + 1]); - } - // Avoid negative target rates comming from too large local reductions. - const Scalar target_rate_available = std::max(Scalar(1e-12), target / efficiencyFactor); - Scalar scale = 1.0; - if (current_rate_available > 1e-12) - scale = target_rate_available / current_rate_available; - - return std::make_pair(current_rate_available > target_rate_available, scale); -} - - -template -Scalar -WellGroupHelpers:: -getWellGroupTargetInjector(const std::string& name, - const std::string& parent, - const Group& group, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - Phase injectionPhase, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - DeferredLogger& deferred_logger) -{ - // This function computes a wells group target. - // 'parent' will be the name of 'group'. But if we recurse, 'name' and - // 'parent' will stay fixed while 'group' will be higher up - // in the group tree. - // Eficiencyfactor is the well efficiency factor for the first group the well is - // part of. Later it is the accumulated factor including the group efficiency factor - // of the child of group. - auto currentGroupControl = group_state.injection_control(group.name(), injectionPhase); - if (currentGroupControl == Group::InjectionCMode::FLD || currentGroupControl == Group::InjectionCMode::NONE) { - // Return if we are not available for parent group. - if (!group.injectionGroupControlAvailable(injectionPhase)) { - return std::numeric_limits::max(); - } - // Otherwise: check production share of parent's control. - const auto& parentGroup = schedule.getGroup(group.parent(), reportStepIdx); - return getWellGroupTargetInjector(name, - parent, - parentGroup, - wellState, - group_state, - reportStepIdx, - guideRate, - rates, - injectionPhase, - efficiencyFactor * group.getGroupEfficiencyFactor(), - schedule, - summaryState, - resv_coeff, - deferred_logger); - } - - // This can be false for FLD-controlled groups, we must therefore - // check for FLD first (done above). - if (!group.isInjectionGroup()) { - return std::numeric_limits::max(); - } - - // If we are here, we are at the topmost group to be visited in the recursion. - // This is the group containing the control we will check against. - Scalar sales_target = 0; - if (schedule[reportStepIdx].gconsale().has(group.name())) { - const auto& gconsale = schedule[reportStepIdx].gconsale().get(group.name(), summaryState); - sales_target = gconsale.sales_target; - } - WGHelpers::InjectionTargetCalculator tcalc(currentGroupControl, - wellState.phaseUsageInfo(), - resv_coeff, - group.name(), - sales_target, - group_state, - injectionPhase, - group.has_gpmaint_control(injectionPhase, - currentGroupControl), - deferred_logger); - - WGHelpers::FractionCalculator fcalc(schedule, - wellState, - group_state, - summaryState, - reportStepIdx, - guideRate, - tcalc.guideTargetMode(), - false, - injectionPhase); - - auto localFraction = [&](const std::string& child, const std::string& always_incluced_name) { return fcalc.localFraction(child, always_incluced_name); }; - - auto localReduction = [&](const std::string& group_name) - { - const std::vector& groupTargetReductions = - group_state.injection_reduction_rates(group_name); - return tcalc.calcModeRateFromRates(groupTargetReductions); - }; - - std::optional ctrl; - if (!group.has_gpmaint_control(injectionPhase, currentGroupControl)) - ctrl = group.injectionControls(injectionPhase, summaryState); - - - const Scalar orig_target = tcalc.groupTarget(ctrl, deferred_logger); - // Assume we have a chain of groups as follows: BOTTOM -> MIDDLE -> TOP. - // Then ... - // TODO finish explanation. - const Scalar current_rate_available - = tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers. - const auto chain = groupChainTopBot(name, group.name(), schedule, reportStepIdx); - // Because 'name' is the last of the elements, and not an ancestor, we subtract one below. - const std::size_t num_ancestors = chain.size() - 1; - // we need to find out the level where the local reduction is applied - std::size_t local_reduction_level = 0; - for (std::size_t ii = 1; ii < num_ancestors; ++ii) { - const int num_gr_ctrl = groupControlledWells(schedule, - wellState, - group_state, - reportStepIdx, - chain[ii], - "", - /*is_producer*/ false, - injectionPhase); - if (guideRate->has(chain[ii], injectionPhase) && num_gr_ctrl > 0) { - local_reduction_level = ii; - } - } - - // Compute portion of target corresponding to current_rate_available - Scalar target = orig_target; - for (std::size_t ii = 0; ii < num_ancestors; ++ii) { - if ((ii == 0) || guideRate->has(chain[ii], injectionPhase)) { - // Apply local reductions only at the control level - // (top) and for levels where we have a specified - // group guide rate. - if (local_reduction_level >= ii ) { - target -= localReduction(chain[ii]); - } - - // If we are under individual control we need to add the wells rate back at the level where it is included in the local reduction - if (local_reduction_level == ii && !wellState.isInjectionGrup(name)) { - target += current_rate_available * efficiencyFactor; - } - } - if (!wellState.isInjectionGrup(name)) { - target *= localFraction(chain[ii + 1], name); - } else { - target *= localFraction(chain[ii + 1], chain[ii + 1]); - } - } - // Avoid negative target rates comming from too large local reductions. - return std::max(Scalar(0.0), target / efficiencyFactor); -} - - -template -std::pair, Scalar> -WellGroupHelpers:: -worstOffendingWell(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const Group::ProductionCMode& offendedControl, - const Parallel::Communication& comm, - const WellState& wellState, - DeferredLogger& deferred_logger) -{ - std::pair, Scalar> offending_well {std::nullopt, 0.0}; - for (const std::string& child_group : group.groups()) { - const auto& this_group = schedule.getGroup(child_group, reportStepIdx); - const auto & offending_well_this = worstOffendingWell(this_group, - schedule, - reportStepIdx, - offendedControl, - comm, - wellState, - deferred_logger); - if (offending_well_this.second > offending_well.second) { - offending_well = offending_well_this; - } - } - - for (const std::string& child_well : group.wells()) { - - const auto& well_index = wellState.index(child_well); - Scalar violating_rate = 0.0; - Scalar prefered_rate = 0.0; - const auto& pu = wellState.phaseUsageInfo(); - if (well_index.has_value() && wellState.wellIsOwned(well_index.value(), child_well)) - { - const auto& ws = wellState.well(child_well); - switch (offendedControl){ - case Group::ProductionCMode::ORAT: { - const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); - violating_rate = ws.surface_rates[oil_pos]; - break; - } - case Group::ProductionCMode::GRAT: { - const int gas_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); - violating_rate = ws.surface_rates[gas_pos]; - break; - } - case Group::ProductionCMode::WRAT: { - const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); - violating_rate = ws.surface_rates[water_pos]; - break; - } - case Group::ProductionCMode::LRAT: { - assert(pu.phaseIsActive(IndexTraits::oilPhaseIdx)); - assert(pu.phaseIsActive(IndexTraits::waterPhaseIdx)); - const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); - const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); - violating_rate = ws.surface_rates[oil_pos] + - ws.surface_rates[water_pos]; - break; - } - case Group::ProductionCMode::RESV: - for (int p = 0; p < wellState.numPhases(); ++p) { - violating_rate += ws.reservoir_rates[p]; - } - break; - case Group::ProductionCMode::NONE: - break; - case Group::ProductionCMode::FLD: - break; - case Group::ProductionCMode::PRBL: - OPM_DEFLOG_THROW(std::runtime_error, - "Group " + group.name() + - "GroupProductionCMode PRBL not implemented", deferred_logger); - break; - case Group::ProductionCMode::CRAT: - OPM_DEFLOG_THROW(std::runtime_error, - "Group " + group.name() + - "GroupProductionCMode CRAT not implemented", deferred_logger); - break; - } - const auto preferred_phase = schedule.getWell(child_well, reportStepIdx).getPreferredPhase(); - switch (preferred_phase) { - case Phase::OIL: { - const int oil_pos = pu.canonicalToActivePhaseIdx(IndexTraits::oilPhaseIdx); - prefered_rate = ws.surface_rates[oil_pos]; - break; - } - case Phase::GAS: { - const int gas_pos = pu.canonicalToActivePhaseIdx(IndexTraits::gasPhaseIdx); - prefered_rate = ws.surface_rates[gas_pos]; - break; - } - case Phase::WATER: { - const int water_pos = pu.canonicalToActivePhaseIdx(IndexTraits::waterPhaseIdx); - prefered_rate = ws.surface_rates[water_pos]; - break; - } - default: - // No others supported. - break; - } - } - violating_rate = comm.sum(violating_rate); - if (violating_rate < 0 ) { // only check producing wells - prefered_rate = comm.sum(prefered_rate); - Scalar fraction = prefered_rate < -1e-16 ? violating_rate / prefered_rate : 1.0; - if ( fraction > offending_well.second) { - offending_well = {child_well, fraction}; - } - } - } - return offending_well; -} - -template -template -void WellGroupHelpers:: -setRegionAveragePressureCalculator(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const FieldPropsManager& fp, - std::map>& regionalAveragePressureCalculator) -{ - for (const std::string& groupName : group.groups()) { - setRegionAveragePressureCalculator( schedule.getGroup(groupName, reportStepIdx), schedule, - reportStepIdx, fp, regionalAveragePressureCalculator); - } - const auto& gpm = group.gpmaint(); - if (!gpm) - return; - - const auto& reg = gpm->region(); - if (!reg) - return; - - if (regionalAveragePressureCalculator.count(reg->first) == 0) { - const std::string name = (reg->first.rfind("FIP", 0) == 0) ? reg->first : "FIP" + reg->first; - const auto& fipnum = fp.get_int(name); - regionalAveragePressureCalculator[reg->first] = std::make_unique(fipnum); - } -} - -#include - -template -using AvgP = RegionAverageCalculator:: - AverageRegionalPressure,std::vector>; - -template -using AvgPMap = std::map>>; - -#define INSTANTIATE_TYPE(T) \ -template class WellGroupHelpers; \ -template void WellGroupHelpers:: \ -updateGpMaintTargetForGroups>(const Group&, \ -const Schedule&, \ -const AvgPMap&, \ -int, \ -double, \ -const WellState&, \ -GroupState&); \ -template void WellGroupHelpers:: \ -setRegionAveragePressureCalculator>(const Group&, \ -const Schedule&, \ -const int, \ -const FieldPropsManager&, \ -AvgPMap&); - -INSTANTIATE_TYPE(double) - -#if FLOW_INSTANTIATE_FLOAT -INSTANTIATE_TYPE(float) -#endif - -} // namespace Opm::WellGroupHelpers diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp deleted file mode 100644 index dce1ad4803..0000000000 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ /dev/null @@ -1,351 +0,0 @@ -/* - Copyright 2019 Norce. - - 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_WELLGROUPHELPERS_HEADER_INCLUDED -#define OPM_WELLGROUPHELPERS_HEADER_INCLUDED - -#include "opm/material/fluidsystems/PhaseUsageInfo.hpp" - -#include -#include -#include -#include -#include - -#include -#include -#include - -namespace Opm { - -class DeferredLogger; -class Group; -template class GroupState; -namespace Network { class ExtNetwork; } -class Schedule; -template class VFPProdProperties; -template class WellState; -class FieldPropsManager; - -namespace Network { class ExtNetwork; } - -template -class WellGroupHelpers -{ -public: - - using WellStateType = WellState; - - static Scalar sumWellPhaseRates(bool res_rates, - const Opm::Group& group, - const Opm::Schedule& schedule, - const WellStateType& wellState, - const SummaryState& summaryState, - const int reportStepIdx, - const int phasePos, - const bool injector, - const bool network = false); - - static Scalar satelliteInjectionRate(const ScheduleState& sched, - const Group& group, - const PhaseUsageInfo& pu, - const int phase_pos, - bool res_rates); - - static Scalar satelliteProductionRate(const SummaryState& summaryState, - const ScheduleState& sched, - const Group& group, - const GSatProd::GSatProdGroupProp::Rate rateComp, - bool res_rates); - - static std::optional - selectRateComponent(const PhaseUsageInfo& pu, const int phasePos); - - static void setCmodeGroup(const Group& group, - const Schedule& schedule, - const SummaryState& summaryState, - const int reportStepIdx, - GroupState& group_state); - - static void accumulateGroupEfficiencyFactor(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - Scalar& factor); - - static Scalar sumWellSurfaceRates(const Group& group, - const Schedule& schedule, - const WellStateType& wellState, - const int reportStepIdx, - const int phasePos, - const bool injector, - const SummaryState& summaryState); - - /// Returns the name of the worst offending well and its fraction (i.e. violated_phase / preferred_phase) - static std::pair, Scalar> - worstOffendingWell(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const Group::ProductionCMode& offendedControl, - const Parallel::Communication& comm, - const WellStateType& wellState, - DeferredLogger& deferred_logger); - - static Scalar sumWellResRates(const Group& group, - const Schedule& schedule, - const WellStateType& wellState, - const int reportStepIdx, - const int phasePos, - const bool injector, - const SummaryState& summaryState); - - static Scalar sumSolventRates(const Group& group, - const Schedule& schedule, - const WellStateType& wellState, - const int reportStepIdx, - const bool injector); - - static void updateGroupTargetReduction(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const bool isInjector, - const GuideRate& guide_rate, - const WellStateType& wellState, - const SummaryState& summaryState, - GroupState& group_state, - std::vector& groupTargetReduction); - - static void updateVREPForGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellStateType& wellState, - GroupState& group_state, - const SummaryState& summaryState); - - template - static void updateGpMaintTargetForGroups(const Group& group, - const Schedule& schedule, - const RegionalValues& regional_values, - const int reportStepIdx, - const double dt, - const WellStateType& well_state, - GroupState& group_state); - - static void updateReservoirRatesInjectionGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellStateType& wellState, - GroupState& group_state, - const SummaryState& summaryState); - - static void updateSurfaceRatesInjectionGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellStateType& wellState, - GroupState& group_state, - const SummaryState& summaryState); - - static void updateWellRates(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellStateType& wellStateNupcol, - WellStateType& wellState); - - static void updateGroupProductionRates(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const WellStateType& wellState, - GroupState& group_state, - const SummaryState& summaryState); - - static void updateNetworkLeafNodeProductionRates(const Schedule& schedule, - const int reportStepIdx, - const WellStateType& wellState, - GroupState& group_state, - const SummaryState& summaryState); - - - static void updateWellRatesFromGroupTargetScale(const Scalar scale, - const Group& group, - const Schedule& schedule, - const int reportStepIdx, - bool isInjector, - const GroupState& group_state, - WellStateType& wellState); - - static void updateREINForGroups(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const SummaryState& st, - const WellStateType& wellState, - GroupState& group_state, - bool sum_rank); - - - static std::map - computeNetworkPressures(const Network::ExtNetwork& network, - const WellStateType& well_state, - const GroupState& group_state, - const VFPProdProperties& vfp_prod_props, - const Schedule& schedule, - const Parallel::Communication& comm, - const int report_time_step); - - static GuideRate::RateVector - getWellRateVector(const WellStateType& well_state, - const std::string& name); - - static GuideRate::RateVector - getProductionGroupRateVector(const GroupState& group_state, - const PhaseUsageInfo& pu, - const std::string& group_name); - - static Scalar getGuideRate(const std::string& name, - const Schedule& schedule, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const GuideRateModel::Target target); - - static Scalar getGuideRateInj(const std::string& name, - const Schedule& schedule, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const GuideRateModel::Target target, - const Phase& injectionPhase); - - /// update the number of wells that are actively under group control for a given group with name given by group_name - /// its main usage is to detect cases where there is no wells under group control - static int updateGroupControlledWells(const Schedule& schedule, - const WellStateType& well_state, - GroupState& group_state, - const SummaryState& summary_state, - const GuideRate* guideRate, - const int report_step, - const std::string& group_name, - const bool is_production_group, - const Phase injection_phase); - - /// returns the number of wells that are actively under group control for a given group with name given by group_name - static int groupControlledWells(const Schedule& schedule, - const WellStateType& well_state, - const GroupState& group_state, - const int report_step, - const std::string& group_name, - const std::string& always_included_child, - const bool is_production_group, - const Phase injection_phase); - - static std::pair - checkGroupConstraintsInj(const std::string& name, - const std::string& parent, - const Group& group, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - Phase injectionPhase, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - const bool check_guide_rate, - DeferredLogger& deferred_logger); - - static Scalar - getWellGroupTargetInjector(const std::string& name, - const std::string& parent, - const Group& group, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - Phase injectionPhase, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - DeferredLogger& deferred_logger); - - static std::vector - groupChainTopBot(const std::string& bottom, - const std::string& top, - const Schedule& schedule, - const int report_step); - - - // check if well/group bottom is a sub well/group of the group top - static bool - isInGroupChainTopBot(const std::string& bottom, - const std::string& top, - const Schedule& schedule, - const int report_step); - - static std::string - control_group(const Group& group, - const GroupState& group_state, - const int reportStepIdx, - const Schedule& schedule); - - static std::pair - checkGroupConstraintsProd(const std::string& name, - const std::string& parent, - const Group& group, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - const bool check_guide_rate, - DeferredLogger& deferred_logger); - static Scalar - getWellGroupTargetProducer(const std::string& name, - const std::string& parent, - const Group& group, - const WellStateType& wellState, - const GroupState& group_state, - const int reportStepIdx, - const GuideRate* guideRate, - const Scalar* rates, - const Scalar efficiencyFactor, - const Schedule& schedule, - const SummaryState& summaryState, - const std::vector& resv_coeff, - DeferredLogger& deferred_logger); - - template - static void setRegionAveragePressureCalculator(const Group& group, - const Schedule& schedule, - const int reportStepIdx, - const FieldPropsManager& fp, - std::map>& regionalAveragePressureCalculator); -}; - -} // namespace Opm - -#endif diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 8e86fc8f26..f3ddb10530 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -69,6 +69,7 @@ namespace Opm class WellInjectionProperties; class WellProductionProperties; +template class WellGroupHelper; template class WellInterface : public WellInterfaceIndices, @@ -98,6 +99,7 @@ public: using WellStateType = WellState; using SingleWellStateType = SingleWellState; + using WellGroupHelperType = WellGroupHelper; using RateConverterType = typename WellInterfaceFluidSystem::RateConverterType; @@ -167,8 +169,8 @@ public: void assembleWellEq(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger); void assembleWellEqWithoutIteration(const Simulator& simulator, @@ -180,8 +182,8 @@ public: // TODO: better name or further refactoring the function to make it more clear void prepareWellBeforeAssembling(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger); virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator, @@ -191,6 +193,7 @@ public: virtual std::optional computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, const Scalar alq_value, DeferredLogger& deferred_logger, @@ -212,11 +215,12 @@ public: // TODO: before we decide to put more information under mutable, this function is not const virtual void computeWellPotentials(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, std::vector& well_potentials, DeferredLogger& deferred_logger) = 0; virtual void updateWellStateWithTarget(const Simulator& simulator, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger) const; @@ -224,6 +228,7 @@ public: virtual void computeWellRatesWithBhpIterations(const Simulator& simulator, const Scalar& bhp, + const WellGroupHelperType& wgHelper, std::vector& well_flux, DeferredLogger& deferred_logger) const = 0; @@ -242,21 +247,22 @@ public: bool updateWellStateWithTHPTargetProd(const Simulator& simulator, WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) const; enum class IndividualOrGroup { Individual, Group, Both }; bool updateWellControl(const Simulator& simulator, const IndividualOrGroup iog, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) /* const */; bool updateWellControlAndStatusLocalIteration(const Simulator& simulator, - WellStateType& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, const Scalar WQTotal, + WellStateType& well_state, DeferredLogger& deferred_logger, const bool fixed_control = false, const bool fixed_status = false); @@ -291,8 +297,8 @@ public: // Simulator is not const is because that assembleWellEq is non-const Simulator void wellTesting(const Simulator& simulator, const double simulation_time, - /* const */ WellStateType& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, + WellStateType& well_state, WellTestState& welltest_state, GLiftEclWells& ecl_well_map, std::map& open_times, @@ -300,6 +306,7 @@ public: void checkWellOperability(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger); void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& simulator, @@ -312,10 +319,11 @@ public: // mostly related to BHP limit and THP limit void updateWellOperability(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger); bool updateWellOperabilityFromWellEq(const Simulator& simulator, - const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger); // update perforation water throughput based on solved water rate @@ -336,8 +344,8 @@ public: DeferredLogger& deferred_logger) const; void solveWellEquation(const Simulator& simulator, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger); const std::vector& connectionRates() const @@ -360,8 +368,8 @@ public: const double dt, const WellInjectionControls& inj_controls, const WellProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger, const bool fixed_control = false, const bool fixed_status = false) = 0; @@ -394,6 +402,7 @@ protected: // check whether the well is operable under THP limit with current reservoir condition virtual void checkOperabilityUnderTHPLimit(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) = 0; virtual void updateIPR(const Simulator& simulator, @@ -412,8 +421,8 @@ protected: const double dt, const WellInjectionControls& inj_controls, const WellProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) = 0; virtual void updateIPRImplicit(const Simulator& simulator, @@ -422,39 +431,42 @@ protected: bool iterateWellEquations(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger); bool solveWellWithOperabilityCheck(const Simulator& simulator, const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger); std::optional estimateOperableBhp(const Simulator& ebos_simulator, const double dt, - WellStateType& well_state, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, + WellStateType& well_state, DeferredLogger& deferred_logger); bool solveWellWithBhp(const Simulator& simulator, const double dt, const Scalar bhp, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger); bool solveWellWithZeroRate(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger); bool solveWellForTesting(const Simulator& simulator, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp index 91878ea45e..6882d6b61a 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp @@ -37,7 +37,7 @@ #include #include #include -#include +#include #include namespace Opm @@ -179,13 +179,16 @@ checkIndividualConstraints(SingleWellState& ws, template bool WellInterfaceFluidSystem:: -checkGroupConstraints(WellState& well_state, - const GroupState& group_state, +checkGroupConstraints(const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, const bool check_guide_rate, + WellStateType& well_state, DeferredLogger& deferred_logger) const { + const auto& group_state = wgHelper.groupState(); + + if (!this->wellEcl().isAvailableForGroupControl()) return false; @@ -201,18 +204,18 @@ checkGroupConstraints(WellState& well_state, this->rateConverter().calcInjCoeff(id, region, coeff); }; - return WellGroupConstraints(*this).checkGroupConstraints(well_state, group_state, + return WellGroupConstraints(*this).checkGroupConstraints(wgHelper, schedule, summaryState, - rCoeff, check_guide_rate, deferred_logger); + rCoeff, check_guide_rate, well_state, deferred_logger); } template bool WellInterfaceFluidSystem:: -checkConstraints(WellState& well_state, - const GroupState& group_state, +checkConstraints(const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, + WellStateType& well_state, DeferredLogger& deferred_logger) const { const bool ind_broken = checkIndividualConstraints(well_state.well(this->index_of_well_), @@ -220,8 +223,8 @@ checkConstraints(WellState& well_state, if (ind_broken) { return true; } else { - return checkGroupConstraints(well_state, group_state, schedule, - summaryState, true, deferred_logger); + return checkGroupConstraints(wgHelper, schedule, + summaryState, true, well_state, deferred_logger); } } diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.hpp b/opm/simulators/wells/WellInterfaceFluidSystem.hpp index 061e26fdcc..2f4654a95c 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.hpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.hpp @@ -43,6 +43,7 @@ class Schedule; struct RatioLimitCheckReport; template class SingleWellState; template class WellState; +template class WellGroupHelper; template class WellInterfaceFluidSystem : public WellInterfaceGeneric @@ -57,6 +58,8 @@ public: using Scalar = typename FluidSystem::Scalar; using IndexTraits = typename FluidSystem::IndexTraitsType; using ModelParameters = typename WellInterfaceGeneric::ModelParameters; + using WellGroupHelperType = WellGroupHelper; + using WellStateType = WellState; static constexpr int Water = IndexTraits::waterPhaseIdx; static constexpr int Oil = IndexTraits::oilPhaseIdx; @@ -88,17 +91,17 @@ protected: const std::optional& inj_controls = std::nullopt, const std::optional& prod_controls = std::nullopt) const; - bool checkGroupConstraints(WellState& well_state, - const GroupState& group_state, + bool checkGroupConstraints(const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, const bool check_guide_rate, + WellStateType& well_state, DeferredLogger& deferred_logger) const; - bool checkConstraints(WellState& well_state, - const GroupState& group_state, + bool checkConstraints(const WellGroupHelperType& wgHelper, const Schedule& schedule, const SummaryState& summaryState, + WellStateType& well_state, DeferredLogger& deferred_logger) const; std::optional diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 06d9493607..f61e178a10 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -188,8 +188,8 @@ namespace Opm WellInterface:: updateWellControl(const Simulator& simulator, const IndividualOrGroup iog, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) /* const */ { OPM_TIMEFUNCTION(); @@ -233,10 +233,12 @@ namespace Opm if (iog == IndividualOrGroup::Individual) { changed = this->checkIndividualConstraints(ws, summaryState, deferred_logger); } else if (iog == IndividualOrGroup::Group) { - changed = this->checkGroupConstraints(well_state, group_state, schedule, summaryState, true, deferred_logger); + changed = this->checkGroupConstraints( + wgHelper, schedule, summaryState, true, well_state, deferred_logger + ); } else { assert(iog == IndividualOrGroup::Both); - changed = this->checkConstraints(well_state, group_state, schedule, summaryState, deferred_logger); + changed = this->checkConstraints(wgHelper, schedule, summaryState, well_state, deferred_logger); } Parallel::Communication cc = simulator.vanguard().grid().comm(); // checking whether control changed @@ -262,7 +264,7 @@ namespace Opm if (iterationIdx >= nupcol || this->well_control_log_.empty()) { this->well_control_log_.push_back(from); } - updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger); + updateWellStateWithTarget(simulator, wgHelper, well_state, deferred_logger); updatePrimaryVariables(simulator, well_state, deferred_logger); } @@ -273,11 +275,11 @@ namespace Opm bool WellInterface:: updateWellControlAndStatusLocalIteration(const Simulator& simulator, - WellStateType& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, const Scalar wqTotal, + WellStateType& well_state, DeferredLogger& deferred_logger, const bool fixed_control, const bool fixed_status) @@ -317,7 +319,9 @@ namespace Opm changed = this->checkIndividualConstraints(ws, summary_state, deferred_logger, inj_controls, prod_controls); } if (hasGroupControl && this->param_.check_group_constraints_inner_well_iterations_) { - changed = changed || this->checkGroupConstraints(well_state, group_state, schedule, summary_state, false, deferred_logger); + changed = changed || this->checkGroupConstraints( + wgHelper, schedule, summary_state, false, well_state, deferred_logger + ); } if (changed) { @@ -327,7 +331,7 @@ namespace Opm ws.thp = this->getTHPConstraint(summary_state); } else { // don't call for thp since this might trigger additional local solve - updateWellStateWithTarget(simulator, group_state, well_state, deferred_logger); + updateWellStateWithTarget(simulator, wgHelper, well_state, deferred_logger); } updatePrimaryVariables(simulator, well_state, deferred_logger); } @@ -382,17 +386,23 @@ namespace Opm WellInterface:: wellTesting(const Simulator& simulator, const double simulation_time, - /* const */ WellStateType& well_state, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, + WellStateType& well_state, WellTestState& well_test_state, GLiftEclWells& ecl_well_map, std::map& open_times, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); + const auto& group_state = wgHelper.groupState(); deferred_logger.info(" well " + this->name() + " is being tested"); + WellGroupHelperType wgHelper_copy = wgHelper; WellStateType well_state_copy = well_state; + // Ensure that wgHelper uses well_state_copy as WellState for the well testing + // and the guard ensures that the original well state is restored at scope exit, i.e. at + // the end of this function. + auto guard = wgHelper_copy.pushWellState(well_state_copy); auto& ws = well_state_copy.well(this->indexOfWell()); const auto& summary_state = simulator.vanguard().summaryState(); @@ -430,7 +440,7 @@ namespace Opm // untill the number of closed completions do not increase anymore. while (testWell) { const std::size_t original_number_closed_completions = welltest_state_temp.num_closed_completions(); - bool converged = solveWellForTesting(simulator, well_state_copy, group_state, deferred_logger); + bool converged = solveWellForTesting(simulator, wgHelper_copy, well_state_copy, deferred_logger); if (!converged) { const auto msg = fmt::format("WTEST: Well {} is not solvable (physical)", this->name()); deferred_logger.debug(msg); @@ -438,7 +448,7 @@ namespace Opm } - updateWellOperability(simulator, well_state_copy, deferred_logger); + updateWellOperability(simulator, well_state_copy, wgHelper_copy, deferred_logger); if ( !this->isOperableAndSolvable() ) { const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name()); deferred_logger.debug(msg); @@ -446,7 +456,7 @@ namespace Opm } std::vector potentials; try { - computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger); + computeWellPotentials(simulator, well_state_copy, wgHelper_copy, potentials, deferred_logger); } catch (const std::exception& e) { const std::string msg = fmt::format("well {}: computeWellPotentials() " "failed during testing for re-opening: ", @@ -502,8 +512,8 @@ namespace Opm WellInterface:: iterateWellEquations(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); @@ -517,12 +527,16 @@ namespace Opm try { // TODO: the following two functions will be refactored to be one to reduce the code duplication if (!this->param_.local_well_solver_control_switching_){ - converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger); } else { if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && (well_state.well(this->index_of_well_).status == WellStatus::OPEN)) { - converged = solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = solveWellWithOperabilityCheck( + simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger + ); } else { - converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = this->iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger + ); } } @@ -567,8 +581,8 @@ namespace Opm const double dt, const Well::InjectionControls& inj_controls, const Well::ProductionControls& prod_controls, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); @@ -578,14 +592,16 @@ namespace Opm // if well is stopped, check if we can reopen if (this->wellIsStopped()) { this->openWell(); - auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger); + auto bhp_target = estimateOperableBhp( + simulator, dt, wgHelper, summary_state, well_state, deferred_logger + ); if (!bhp_target.has_value()) { // no intersection with ipr const auto msg = fmt::format("estimateOperableBhp: Did not find operable BHP for well {}", this->name()); deferred_logger.debug(msg); // well can't operate using explicit fractions stop the well // solve with zero rates - converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger); + converged = solveWellWithZeroRate(simulator, dt, wgHelper, well_state, deferred_logger); this->stopWell(); this->operability_status_.can_obtain_bhp_with_thp_limit = false; this->operability_status_.obey_thp_limit_under_bhp_limit = false; @@ -595,11 +611,13 @@ namespace Opm ws.thp = this->getTHPConstraint(summary_state); const Scalar bhp = std::max(bhp_target.value(), static_cast(prod_controls.bhp_limit)); - solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger); + solveWellWithBhp(simulator, dt, bhp, wgHelper, well_state, deferred_logger); } } // solve well-equation - converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = this->iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger + ); const bool isThp = ws.production_cmode == Well::ProducerCMode::THP; @@ -619,10 +637,14 @@ namespace Opm if (bhp_stable.has_value() && cur_bhp - bhp_stable.value() > cur_bhp*reltol){ const auto msg = fmt::format("Well {} converged to an unstable solution, re-solving", this->name()); deferred_logger.debug(msg); - solveWellWithBhp(simulator, dt, bhp_stable.value(), well_state, deferred_logger); + solveWellWithBhp( + simulator, dt, bhp_stable.value(), wgHelper, well_state, deferred_logger + ); // re-solve with hopefully good initial guess ws.thp = this->getTHPConstraint(summary_state); - converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = this->iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger + ); } } } @@ -631,11 +653,13 @@ namespace Opm // Well did not converge, switch to explicit fractions this->operability_status_.use_vfpexplicit = true; this->openWell(); - auto bhp_target = estimateOperableBhp(simulator, dt, well_state, summary_state, deferred_logger); + auto bhp_target = estimateOperableBhp( + simulator, dt, wgHelper, summary_state, well_state, deferred_logger + ); if (!bhp_target.has_value()) { // solve with zero rate // well can't operate using explicit fractions stop the well - converged = solveWellWithZeroRate(simulator, dt, well_state, deferred_logger); + converged = solveWellWithZeroRate(simulator, dt, wgHelper, well_state, deferred_logger); this->stopWell(); this->operability_status_.can_obtain_bhp_with_thp_limit = false; this->operability_status_.obey_thp_limit_under_bhp_limit = false; @@ -644,15 +668,17 @@ namespace Opm // solve well with the estimated target bhp (or limit) const Scalar bhp = std::max(bhp_target.value(), static_cast(prod_controls.bhp_limit)); - solveWellWithBhp(simulator, dt, bhp, well_state, deferred_logger); + solveWellWithBhp( + simulator, dt, bhp, wgHelper, well_state, deferred_logger + ); ws.thp = this->getTHPConstraint(summary_state); const auto msg = fmt::format("Well {} did not converge, re-solving with explicit fractions for VFP caculations.", this->name()); deferred_logger.debug(msg); converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, + wgHelper, well_state, - group_state, deferred_logger); } } @@ -667,13 +693,16 @@ namespace Opm WellInterface:: estimateOperableBhp(const Simulator& simulator, const double dt, - WellStateType& well_state, + const WellGroupHelperType& wgHelper, const SummaryState& summary_state, + WellStateType& well_state, DeferredLogger& deferred_logger) { if (!this->wellHasTHPConstraints(summary_state)) { const Scalar bhp_limit = WellBhpThpCalculator(*this).mostStrictBhpFromBhpLimits(summary_state); - const bool converged = solveWellWithBhp(simulator, dt, bhp_limit, well_state, deferred_logger); + const bool converged = solveWellWithBhp( + simulator, dt, bhp_limit, wgHelper, well_state, deferred_logger + ); if (!converged || this->wellIsStopped()) { return std::nullopt; } @@ -685,7 +714,9 @@ namespace Opm // Get minimal bhp from vfp-curve Scalar bhp_min = WellBhpThpCalculator(*this).calculateMinimumBhpFromThp(well_state, this->well_ecl_, summary_state, this->getRefDensity()); // Solve - const bool converged = solveWellWithBhp(simulator, dt, bhp_min, well_state, deferred_logger); + const bool converged = solveWellWithBhp( + simulator, dt, bhp_min, wgHelper, well_state, deferred_logger + ); if (!converged || this->wellIsStopped()) { return std::nullopt; } @@ -701,12 +732,19 @@ namespace Opm solveWellWithBhp(const Simulator& simulator, const double dt, const Scalar bhp, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); // Solve a well using single bhp-constraint (but close if not operable under this) auto group_state = GroupState(); // empty group + WellGroupHelperType wgHelper_copy = wgHelper; + // Ensure that wgHelper_copy uses the empty group state as GroupState for iterateWellEqWithSwitching() + // and the guard ensures that the original group state is restored at scope exit, i.e. at + // the end of this function. + auto group_guard = wgHelper_copy.pushGroupState(group_state); + auto inj_controls = Well::InjectionControls(0); auto prod_controls = Well::ProductionControls(0); auto& ws = well_state.well(this->index_of_well_); @@ -726,7 +764,10 @@ namespace Opm // update well-state ws.bhp = bhp; // solve - const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger, /*fixed_control*/true); + const bool converged = this->iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper_copy, + well_state, deferred_logger, /*fixed_control*/true + ); ws.injection_cmode = cmode_inj; ws.production_cmode = cmode_prod; return converged; @@ -737,6 +778,7 @@ namespace Opm WellInterface:: solveWellWithZeroRate(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger) { @@ -746,9 +788,18 @@ namespace Opm this->stopWell(); auto group_state = GroupState(); // empty group + WellGroupHelperType wgHelper_copy = wgHelper; + // Ensure that wgHelper_copy uses the empty group state as GroupState for iterateWellEqWithSwitching() + // and the guard ensures that the original group state is restored at scope exit, i.e. at + // the end of this function. + auto group_guard = wgHelper_copy.pushGroupState(group_state); + auto inj_controls = Well::InjectionControls(0); auto prod_controls = Well::ProductionControls(0); - const bool converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger, /*fixed_control*/true, /*fixed_status*/ true); + const bool converged = this->iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper_copy, well_state, + deferred_logger, /*fixed_control*/true, /*fixed_status*/ true + ); this->wellStatus_ = well_status_orig; return converged; } @@ -757,8 +808,8 @@ namespace Opm bool WellInterface:: solveWellForTesting(const Simulator& simulator, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); @@ -773,12 +824,18 @@ namespace Opm try { // TODO: the following two functions will be refactored to be one to reduce the code duplication if (!this->param_.local_well_solver_control_switching_){ - converged = this->iterateWellEqWithControl(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = this->iterateWellEqWithControl( + simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger + ); } else { if (this->param_.use_implicit_ipr_ && this->well_ecl_.isProducer() && (well_state.well(this->index_of_well_).status == WellStatus::OPEN)) { - converged = solveWellWithOperabilityCheck(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = this->solveWellWithOperabilityCheck( + simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger + ); } else { - converged = this->iterateWellEqWithSwitching(simulator, dt, inj_controls, prod_controls, well_state, group_state, deferred_logger); + converged = this->iterateWellEqWithSwitching( + simulator, dt, inj_controls, prod_controls, wgHelper, well_state, deferred_logger + ); } } @@ -803,8 +860,8 @@ namespace Opm void WellInterface:: solveWellEquation(const Simulator& simulator, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); @@ -814,7 +871,7 @@ namespace Opm // keep a copy of the original well state const WellStateType well_state0 = well_state; const double dt = simulator.timeStepSize(); - bool converged = iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger); + bool converged = iterateWellEquations(simulator, dt, wgHelper, well_state, deferred_logger); // Newly opened wells with THP control sometimes struggles to // converge due to bad initial guess. Or due to the simple fact @@ -847,7 +904,7 @@ namespace Opm const std::string msg = std::string("The newly opened well ") + this->name() + std::string(" with THP control did not converge during inner iterations, we try again with bhp control"); deferred_logger.debug(msg); - converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger); + converged = this->iterateWellEquations(simulator, dt, wgHelper, well_state, deferred_logger); } } @@ -866,12 +923,13 @@ namespace Opm WellInterface:: assembleWellEq(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); - prepareWellBeforeAssembling(simulator, dt, well_state, group_state, deferred_logger); + const auto& group_state = wgHelper.groupState(); + prepareWellBeforeAssembling(simulator, dt, wgHelper, well_state, deferred_logger); assembleWellEqWithoutIteration(simulator, dt, well_state, group_state, deferred_logger); } @@ -902,15 +960,15 @@ namespace Opm WellInterface:: prepareWellBeforeAssembling(const Simulator& simulator, const double dt, + const WellGroupHelperType& wgHelper, WellStateType& well_state, - const GroupState& group_state, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); const bool old_well_operable = this->operability_status_.isOperableAndSolvable(); if (this->param_.check_well_operability_iter_) - checkWellOperability(simulator, well_state, deferred_logger); + checkWellOperability(simulator, well_state, wgHelper, deferred_logger); // only use inner well iterations for the first newton iterations. const int iteration_idx = simulator.model().newtonMethod().numIterations(); @@ -933,7 +991,9 @@ namespace Opm } this->stopWell(); changed_to_stopped_this_step_ = true; - bool converged_zero_rate = this->solveWellWithZeroRate(simulator, dt, well_state, deferred_logger); + bool converged_zero_rate = this->solveWellWithZeroRate( + simulator, dt, wgHelper, well_state, deferred_logger + ); if (this->param_.shut_unsolvable_wells_ && !converged_zero_rate ) { this->operability_status_.solvable = false; } @@ -941,7 +1001,9 @@ namespace Opm number_of_well_reopenings_++; return; } - bool converged = this->iterateWellEquations(simulator, dt, well_state, group_state, deferred_logger); + bool converged = this->iterateWellEquations( + simulator, dt, wgHelper, well_state, deferred_logger + ); if (converged) { const bool zero_target = this->wellUnderZeroRateTarget(simulator, well_state, deferred_logger); @@ -965,7 +1027,7 @@ namespace Opm auto well_state_copy = well_state; std::vector potentials; try { - computeWellPotentials(simulator, well_state_copy, potentials, deferred_logger); + computeWellPotentials(simulator, well_state_copy, wgHelper, potentials, deferred_logger); } catch (const std::exception& e) { const std::string msg = fmt::format("well {}: computeWellPotentials() failed " "during attempt to recompute potentials for well: ", @@ -1038,6 +1100,7 @@ namespace Opm WellInterface:: checkWellOperability(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); @@ -1049,12 +1112,12 @@ namespace Opm return; } - updateWellOperability(simulator, well_state, deferred_logger); + updateWellOperability(simulator, well_state, wgHelper, deferred_logger); if (!this->operability_status_.isOperableAndSolvable()) { this->operability_status_.use_vfpexplicit = true; deferred_logger.debug("EXPLICIT_LOOKUP_VFP", "well not operable, trying with explicit vfp lookup: " + this->name()); - updateWellOperability(simulator, well_state, deferred_logger); + updateWellOperability(simulator, well_state, wgHelper, deferred_logger); } } @@ -1127,11 +1190,12 @@ namespace Opm WellInterface:: updateWellOperability(const Simulator& simulator, const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); if (this->param_.local_well_solver_control_switching_) { - const bool success = updateWellOperabilityFromWellEq(simulator, well_state, deferred_logger); + const bool success = updateWellOperabilityFromWellEq(simulator, wgHelper, deferred_logger); if (!success) { this->operability_status_.solvable = false; deferred_logger.debug("Operability check using well equations did not converge for well " @@ -1155,7 +1219,7 @@ namespace Opm } // we do some extra checking for wells under THP control. if (check_thp) { - checkOperabilityUnderTHPLimit(simulator, well_state, deferred_logger); + checkOperabilityUnderTHPLimit(simulator, well_state, wgHelper, deferred_logger); } } @@ -1163,18 +1227,22 @@ namespace Opm bool WellInterface:: updateWellOperabilityFromWellEq(const Simulator& simulator, - const WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) { OPM_TIMEFUNCTION(); // only makes sense if we're using this parameter is true assert(this->param_.local_well_solver_control_switching_); this->operability_status_.resetOperability(); - WellStateType well_state_copy = well_state; - const auto& group_state = simulator.problem().wellModel().groupState(); + WellGroupHelperType wgHelper_copy = wgHelper; + WellStateType well_state_copy = wgHelper_copy.wellState(); const double dt = simulator.timeStepSize(); + // Ensure that wgHelper uses well_state_copy as WellState for iterateWellEquations() + // and the guard ensures that the original well state is restored at scope exit, i.e. at + // the end of this function. + auto guard = wgHelper_copy.pushWellState(well_state_copy); // equations should be converged at this stage, so only one it is needed - bool converged = iterateWellEquations(simulator, dt, well_state_copy, group_state, deferred_logger); + bool converged = iterateWellEquations(simulator, dt, wgHelper_copy, well_state_copy, deferred_logger); return converged; } @@ -1190,11 +1258,12 @@ namespace Opm void WellInterface:: updateWellStateWithTarget(const Simulator& simulator, - const GroupState& group_state, + const WellGroupHelperType& wgHelper, WellStateType& well_state, DeferredLogger& deferred_logger) const { OPM_TIMEFUNCTION(); + const auto& group_state = wgHelper.groupState(); // only bhp and wellRates are used to initilize the primaryvariables for standard wells const auto& well = this->well_ecl_; const int well_index = this->index_of_well_; @@ -1515,7 +1584,7 @@ namespace Opm } case Well::ProducerCMode::THP: { - const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, deferred_logger); + const bool update_success = updateWellStateWithTHPTargetProd(simulator, well_state, wgHelper, deferred_logger); if (!update_success) { // the following is the original way of initializing well state with THP constraint @@ -2037,18 +2106,19 @@ namespace Opm WellInterface:: updateWellStateWithTHPTargetProd(const Simulator& simulator, WellStateType& well_state, + const WellGroupHelperType& wgHelper, DeferredLogger& deferred_logger) const { OPM_TIMEFUNCTION(); const auto& summary_state = simulator.vanguard().summaryState(); auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq( - simulator, summary_state, this->getALQ(well_state), deferred_logger, /*iterate_if_no_solution */ false); + simulator, wgHelper, summary_state, this->getALQ(well_state), deferred_logger, /*iterate_if_no_solution */ false); if (bhp_at_thp_limit) { std::vector rates(this->number_of_phases_, 0.0); if (thp_update_iterations) { computeWellRatesWithBhpIterations(simulator, *bhp_at_thp_limit, - rates, deferred_logger); + wgHelper, rates, deferred_logger); } else { computeWellRatesWithBhp(simulator, *bhp_at_thp_limit, rates, deferred_logger);