diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 4ad30d701..b0da7681e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -309,6 +309,8 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/wells/PerfData.hpp opm/simulators/wells/PerforationData.hpp opm/simulators/wells/RateConverter.hpp + opm/simulators/wells/RegionAttributeHelpers.hpp + opm/simulators/wells/RegionAverageCalculator.hpp opm/simulators/utils/readDeck.hpp opm/simulators/wells/SingleWellState.hpp opm/simulators/wells/TargetCalculator.hpp diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index f1e315827..7399c1de8 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -60,6 +60,7 @@ #include #include #include +#include #include #include #include @@ -132,6 +133,10 @@ namespace Opm { using RateConverterType = RateConverter:: SurfaceToReservoirVoidage >; + // For computing average pressured used by gpmaint + using AverageRegionalPressureType = RegionAverageCalculator:: + AverageRegionalPressure >; + BlackoilWellModel(Simulator& ebosSimulator); void init(); @@ -317,6 +322,8 @@ namespace Opm { bool alternative_well_rate_init_{}; std::unique_ptr rateConverter_{}; + std::unique_ptr regionalAveragePressureCalculator_{}; + SimulatorReportSingle last_report_{}; diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 59af13092..c5c34a49e 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -450,10 +450,15 @@ checkGroupInjectionConstraints(const Group& group, current_rate = comm_.sum(current_rate); const auto& controls = group.injectionControls(phase, this->summaryState_); - if (controls.surface_max_rate < current_rate) { + double target = controls.surface_max_rate; + + if (group.has_gpmaint_control(phase, Group::InjectionCMode::RATE) && this->groupState().has_gpmaint_target(group.name())) + target = this->groupState().gpmaint_target(group.name()); + + if (target < current_rate) { double scale = 1.0; if (current_rate > 1e-12) - scale = controls.surface_max_rate / current_rate; + scale = target / current_rate; return std::make_pair(Group::InjectionCMode::RATE, scale); } } @@ -468,10 +473,15 @@ checkGroupInjectionConstraints(const Group& group, current_rate = comm_.sum(current_rate); const auto& controls = group.injectionControls(phase, this->summaryState_); - if (controls.resv_max_rate < current_rate) { + double target = controls.resv_max_rate; + + if (group.has_gpmaint_control(phase, Group::InjectionCMode::RESV) && this->groupState().has_gpmaint_target(group.name())) + target = this->groupState().gpmaint_target(group.name()); + + if (target < current_rate) { double scale = 1.0; if (current_rate > 1e-12) - scale = controls.resv_max_rate / current_rate; + scale = target / current_rate; return std::make_pair(Group::InjectionCMode::RESV, scale); } } @@ -637,10 +647,14 @@ checkGroupProductionConstraints(const Group& group, // sum over all nodes current_rate = comm_.sum(current_rate); - if (controls.resv_target < current_rate ) { + double target = controls.resv_target; + if (group.has_gpmaint_control(Group::ProductionCMode::RESV) && this->groupState().has_gpmaint_target(group.name())) + target = this->groupState().gpmaint_target(group.name()); + + if ( target < current_rate ) { double scale = 1.0; if (current_rate > 1e-12) - scale = controls.resv_target / current_rate; + scale = target / current_rate; return std::make_pair(Group::ProductionCMode::RESV, scale); } } @@ -1521,6 +1535,8 @@ updateAndCommunicateGroupData(const int reportStepIdx, WellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); WellGroupHelpers::updateReservoirRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); + WellGroupHelpers::updateSurfaceRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, this->groupState()); + WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState()); // We use the rates from the previous time-step to reduce oscillations diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index aa86cc1e2..e016f715e 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -203,10 +203,15 @@ namespace Opm { WellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, this->wellState(), this->groupState()); // Compute reservoir volumes for RESV controls. - rateConverter_.reset(new RateConverterType (phase_usage_, - std::vector(local_num_cells_, 0))); + rateConverter_.reset(new RateConverterType (phase_usage_, std::vector(local_num_cells_, 0))); rateConverter_->template defineState(ebosSimulator_); + // Compute regional average pressures used by gpmaint + const auto& fp = this->eclState_.fieldProps(); + const auto& fipnum = fp.get_int("FIPNUM"); + regionalAveragePressureCalculator_.reset(new AverageRegionalPressureType (phase_usage_,fipnum)); + regionalAveragePressureCalculator_->template defineState(ebosSimulator_); + { const auto& sched_state = this->schedule()[timeStepIdx]; // update VFP properties @@ -458,6 +463,12 @@ namespace Opm { // update the rate converter with current averages pressures etc in rateConverter_->template defineState(ebosSimulator_); + regionalAveragePressureCalculator_->template defineState(ebosSimulator_); + const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx); + WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup, + schedule_, *regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState()); + + // calculate the well potentials try { updateWellPotentials(reportStepIdx, @@ -472,7 +483,6 @@ namespace Opm { updateWellTestState(simulationTime, wellTestState_); // check group sales limits at the end of the timestep - const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx); checkGconsaleLimits(fieldGroup, this->wellState(), ebosSimulator_.episodeIndex(), local_deferredLogger); diff --git a/opm/simulators/wells/GroupState.cpp b/opm/simulators/wells/GroupState.cpp index f663f5e99..ae7be4672 100644 --- a/opm/simulators/wells/GroupState.cpp +++ b/opm/simulators/wells/GroupState.cpp @@ -107,6 +107,27 @@ const std::vector& GroupState::injection_reduction_rates(const std::stri return group_iter->second; } +//------------------------------------------------------------------------- + +bool GroupState::has_injection_surface_rates(const std::string& gname) const { + auto group_iter = this->inj_surface_rates.find(gname); + return (group_iter != this->inj_surface_rates.end()); +} + +void GroupState::update_injection_surface_rates(const std::string& gname, const std::vector& rates) { + if (rates.size() != this->num_phases) + throw std::logic_error("Wrong number of phases"); + + this->inj_surface_rates[gname] = rates; +} + +const std::vector& GroupState::injection_surface_rates(const std::string& gname) const { + auto group_iter = this->inj_surface_rates.find(gname); + if (group_iter == this->inj_surface_rates.end()) + throw std::logic_error("No such group"); + + return group_iter->second; +} //------------------------------------------------------------------------- @@ -246,6 +267,24 @@ GPMaint::State& GroupState::gpmaint(const std::string& gname) { } +//------------------------------------------------------------------------- + +void GroupState::update_gpmaint_target(const std::string& gname, double target) { + this->m_gpmaint_target[gname] = target; +} + +double GroupState::gpmaint_target(const std::string& gname) const { + auto group_iter = this->m_gpmaint_target.find(gname); + if (group_iter == this->m_gpmaint_target.end()) + throw std::logic_error("No such group"); + + return group_iter->second; +} + +bool GroupState::has_gpmaint_target(const std::string& gname) const { + return (this->m_gpmaint_target.count(gname) > 0); +} + //------------------------------------------------------------------------- namespace { diff --git a/opm/simulators/wells/GroupState.hpp b/opm/simulators/wells/GroupState.hpp index e8504963b..2c53e618b 100644 --- a/opm/simulators/wells/GroupState.hpp +++ b/opm/simulators/wells/GroupState.hpp @@ -52,6 +52,11 @@ public: void update_injection_reservoir_rates(const std::string& gname, const std::vector& rates); const std::vector& injection_reservoir_rates(const std::string& gname) const; + bool has_injection_surface_rates(const std::string& gname) const; + void update_injection_surface_rates(const std::string& gname, const std::vector& rates); + const std::vector& injection_surface_rates(const std::string& gname) const; + + void update_injection_rein_rates(const std::string& gname, const std::vector& rates); const std::vector& injection_rein_rates(const std::string& gname) const; @@ -65,6 +70,10 @@ public: double grat_sales_target(const std::string& gname) const; bool has_grat_sales_target(const std::string& gname) const; + void update_gpmaint_target(const std::string& gname, double target); + double gpmaint_target(const std::string& gname) const; + bool has_gpmaint_target(const std::string& gname) const; + bool has_production_control(const std::string& gname) const; void production_control(const std::string& gname, Group::ProductionCMode cmode); Group::ProductionCMode production_control(const std::string& gname) const; @@ -161,11 +170,13 @@ private: std::map production_controls; std::map> prod_red_rates; std::map> inj_red_rates; + std::map> inj_surface_rates; std::map> inj_resv_rates; std::map> inj_potentials; std::map> inj_rein_rates; std::map inj_vrep_rate; std::map m_grat_sales_target; + std::map m_gpmaint_target; std::map, Group::InjectionCMode> injection_controls; diff --git a/opm/simulators/wells/RateConverter.hpp b/opm/simulators/wells/RateConverter.hpp index 01de81bcd..2664d81ed 100644 --- a/opm/simulators/wells/RateConverter.hpp +++ b/opm/simulators/wells/RateConverter.hpp @@ -25,7 +25,7 @@ #include #include #include - +#include #include #include #include @@ -47,345 +47,6 @@ namespace Opm { namespace RateConverter { - /** - * Convenience tools for implementing the rate conversion - * facility. - */ - namespace Details { - namespace Select { - template - struct RegionIDParameter - { - using type = - typename std::remove_reference::type &; - }; - - template - struct RegionIDParameter - { - using type = RegionID; - }; - } // Select - - /** - * \brief Computes the temperature, pressure, and counter increment. - * - * In a parallel run only cells owned contribute to the cell average. - * \tparam is_parallel Whether this is a parallel run. - */ - template - struct AverageIncrementCalculator - { - /** - * \brief Computes the temperature, pressure, and counter increment. - * \param pressure The pressure. - * \param temperature The temperature. - * \param rs The rs. - * \param rv The rv. - * \param cell The current cell index. - * \param ownership A vector indicating whether a cell is owned - * by this process (value 1), or not (value 0). - * \param cell The cell index. - */ - std::tuple - operator()(const std::vector& pressure, - const std::vector& temperature, - const std::vector& rs, - const std::vector& rv, - const std::vector& ownership, - std::size_t cell){ - if ( ownership[cell] ) - { - return std::make_tuple(pressure[cell], - temperature[cell], - rs[cell], - rv[cell], - 1); - } - else - { - return std::make_tuple(0, 0, 0, 0, 0); - } - } - }; - template<> - struct AverageIncrementCalculator - { - std::tuple - operator()(const std::vector& pressure, - const std::vector& temperature, - const std::vector& rs, - const std::vector& rv, - const std::vector&, - std::size_t cell){ - return std::make_tuple(pressure[cell], - temperature[cell], - rs[cell], - rv[cell], - 1); - } - }; - /** - * Provide mapping from Region IDs to user-specified collection - * of per-region attributes. - * - * \tparam RegionId Region identifier type. Must be hashable by - * \code std::hash<> \endcode. Typically a built-in integer - * type--e.g., \c int. - * - * \tparam Attributes User-defined type that represents - * collection of attributes that have meaning in a per-region - * aggregate sense. Must be copy-constructible. - */ - template - class RegionAttributes - { - public: - /** - * Expose \c RegionId as a vocabulary type for use in query - * methods. - */ - using RegionID = - typename Select::RegionIDParameter - ::value>::type; - - /** - * Constructor. - * - * \tparam RMap Class type that implements the RegionMapping - * protocol. Typically an instantiation of \code - * Opm::RegionMapping<> \endcode. - * - * \param[in] rmap Specific region mapping that provides - * reverse lookup from regions to cells. - * - * \param[in] attr Pre-constructed initialiser for \c - * Attributes. - */ - template - RegionAttributes(const RMap& rmap, - const Attributes& attr) - { - using VT = typename AttributeMap::value_type; - - for (const auto& r : rmap.activeRegions()) { - auto v = std::make_unique(attr); - - const auto stat = attr_.insert(VT(r, std::move(v))); - - if (stat.second) { - // New value inserted. - const auto& cells = rmap.cells(r); - - assert (! cells.empty()); - - // Region's representative cell. - stat.first->second->cell_ = cells[0]; - } - } - } - - /** - * Retrieve representative cell in region. - * - * \param[in] reg Specific region. - * - * \return Representative cell in region \p reg. - */ - int cell(const RegionID reg) const - { - return this->find(reg).cell_; - } - - /** - * Request read-only access to region's attributes. - * - * \param[in] reg Specific region. - * - * \return Read-only access to region \p reg's per-region - * attributes. - */ - const Attributes& attributes(const RegionID reg) const - { - return this->find(reg).attr_; - } - - /** - * Request modifiable access to region's attributes. - * - * \param[in] reg Specific region. - * - * \return Read-write access to region \p reg's per-region - * attributes. - */ - Attributes& attributes(const RegionID reg) - { - return this->find(reg).attr_; - } - - private: - /** - * Aggregate per-region attributes along with region's - * representative cell. - */ - struct Value { - Value(const Attributes& attr) - : attr_(attr) - , cell_(-1) - {} - - Attributes attr_; - int cell_; - }; - - using ID = - typename std::remove_reference::type; - - using AttributeMap = - std::unordered_map>; - - AttributeMap attr_; - - /** - * Read-only access to region's properties. - */ - const Value& find(const RegionID reg) const - { - const auto& i = attr_.find(reg); - - if (i == attr_.end()) { - throw std::invalid_argument("Unknown region ID"); - } - - return *i->second; - } - - /** - * Read-write access to region's properties. - */ - Value& find(const RegionID reg) - { - const auto& i = attr_.find(reg); - - if (i == attr_.end()) { - throw std::invalid_argument("Unknown region ID"); - } - - return *i->second; - } - }; - - /** - * Convenience functions for querying presence/absence of - * active phases. - */ - namespace PhaseUsed { - /** - * Active water predicate. - * - * \param[in] pu Active phase object. - * - * \return Whether or not water is an active phase. - */ - inline bool - water(const PhaseUsage& pu) - { - return pu.phase_used[ BlackoilPhases::Aqua ] != 0; - } - - /** - * Active oil predicate. - * - * \param[in] pu Active phase object. - * - * \return Whether or not oil is an active phase. - */ - inline bool - oil(const PhaseUsage& pu) - { - return pu.phase_used[ BlackoilPhases::Liquid ] != 0; - } - - /** - * Active gas predicate. - * - * \param[in] pu Active phase object. - * - * \return Whether or not gas is an active phase. - */ - inline bool - gas(const PhaseUsage& pu) - { - return pu.phase_used[ BlackoilPhases::Vapour ] != 0; - } - } // namespace PhaseUsed - - /** - * Convenience functions for querying numerical IDs - * ("positions") of active phases. - */ - namespace PhasePos { - /** - * Numerical ID of active water phase. - * - * \param[in] pu Active phase object. - * - * \return Non-negative index/position of water if - * active, -1 if not. - */ - inline int - water(const PhaseUsage& pu) - { - int p = -1; - - if (PhaseUsed::water(pu)) { - p = pu.phase_pos[ BlackoilPhases::Aqua ]; - } - - return p; - } - - /** - * Numerical ID of active oil phase. - * - * \param[in] pu Active phase object. - * - * \return Non-negative index/position of oil if - * active, -1 if not. - */ - inline int - oil(const PhaseUsage& pu) - { - int p = -1; - - if (PhaseUsed::oil(pu)) { - p = pu.phase_pos[ BlackoilPhases::Liquid ]; - } - - return p; - } - - /** - * Numerical ID of active gas phase. - * - * \param[in] pu Active phase object. - * - * \return Non-negative index/position of gas if - * active, -1 if not. - */ - inline int - gas(const PhaseUsage& pu) - { - int p = -1; - - if (PhaseUsed::gas(pu)) { - p = pu.phase_pos[ BlackoilPhases::Vapour ]; - } - - return p; - } - } // namespace PhasePos - } // namespace Details /** * Convert component rates at surface conditions to phase @@ -484,7 +145,7 @@ namespace Opm { // only count oil and gas filled parts of the domain double hydrocarbon = 1.0; const auto& pu = phaseUsage_; - if (Details::PhaseUsed::water(pu)) { + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value(); } @@ -496,15 +157,15 @@ namespace Opm { if (hydrocarbonPV > 0.) { auto& attr = attributes_hpv[reg]; attr.pv += hydrocarbonPV; - if (Details::PhaseUsed::oil(pu) && Details::PhaseUsed::gas(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu) && RegionAttributeHelpers::PhaseUsed::gas(pu)) { attr.rs += fs.Rs().value() * hydrocarbonPV; attr.rv += fs.Rv().value() * hydrocarbonPV; } - if (Details::PhaseUsed::oil(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV; attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV; } else { - assert(Details::PhaseUsed::gas(pu)); + assert(RegionAttributeHelpers::PhaseUsed::gas(pu)); attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV; attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV; } @@ -514,18 +175,18 @@ namespace Opm { if (pv_cell > 0.) { auto& attr = attributes_pv[reg]; attr.pv += pv_cell; - if (Details::PhaseUsed::oil(pu) && Details::PhaseUsed::gas(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu) && RegionAttributeHelpers::PhaseUsed::gas(pu)) { attr.rs += fs.Rs().value() * pv_cell; attr.rv += fs.Rv().value() * pv_cell; } - if (Details::PhaseUsed::oil(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell; attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * pv_cell; - } else if (Details::PhaseUsed::gas(pu)) { + } else if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell; attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * pv_cell; } else { - assert(Details::PhaseUsed::water(pu)); + assert(RegionAttributeHelpers::PhaseUsed::water(pu)); attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell; attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell; } @@ -616,13 +277,13 @@ namespace Opm { const double T = ra.temperature; const double saltConcentration = ra.saltConcentration; - const int iw = Details::PhasePos::water(pu); - const int io = Details::PhasePos::oil (pu); - const int ig = Details::PhasePos::gas (pu); + const int iw = RegionAttributeHelpers::PhasePos::water(pu); + const int io = RegionAttributeHelpers::PhasePos::oil (pu); + const int ig = RegionAttributeHelpers::PhasePos::gas (pu); std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0); - if (Details::PhaseUsed::water(pu)) { + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { // q[w]_r = q[w]_s / bw const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration); @@ -637,7 +298,7 @@ namespace Opm { // Determinant of 'R' matrix const double detR = 1.0 - (Rs * Rv); - if (Details::PhaseUsed::oil(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); @@ -645,12 +306,12 @@ namespace Opm { coeff[io] += 1.0 / den; - if (Details::PhaseUsed::gas(pu)) { + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { coeff[ig] -= ra.rv / den; } } - if (Details::PhaseUsed::gas(pu)) { + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv); @@ -658,7 +319,7 @@ namespace Opm { coeff[ig] += 1.0 / den; - if (Details::PhaseUsed::oil(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { coeff[io] -= ra.rs / den; } } @@ -674,13 +335,13 @@ namespace Opm { const double T = ra.temperature; const double saltConcentration = ra.saltConcentration; - const int iw = Details::PhasePos::water(pu); - const int io = Details::PhasePos::oil (pu); - const int ig = Details::PhasePos::gas (pu); + const int iw = RegionAttributeHelpers::PhasePos::water(pu); + const int io = RegionAttributeHelpers::PhasePos::oil (pu); + const int ig = RegionAttributeHelpers::PhasePos::gas (pu); std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0); - if (Details::PhaseUsed::water(pu)) { + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { // q[w]_r = q[w]_s / bw const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration); @@ -688,12 +349,12 @@ namespace Opm { coeff[iw] = 1.0 / bw; } - if (Details::PhaseUsed::oil(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0); coeff[io] += 1.0 / bo; } - if (Details::PhaseUsed::gas(pu)) { + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0); coeff[ig] += 1.0 / bg; } @@ -732,11 +393,11 @@ namespace Opm { const double T = ra.temperature; const double saltConcentration = ra.saltConcentration; - const int iw = Details::PhasePos::water(pu); - const int io = Details::PhasePos::oil (pu); - const int ig = Details::PhasePos::gas (pu); + const int iw = RegionAttributeHelpers::PhasePos::water(pu); + const int io = RegionAttributeHelpers::PhasePos::oil (pu); + const int ig = RegionAttributeHelpers::PhasePos::gas (pu); - if (Details::PhaseUsed::water(pu)) { + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { // q[w]_r = q[w]_s / bw const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration); @@ -764,7 +425,7 @@ namespace Opm { // Determinant of 'R' matrix const double detR = 1.0 - (Rs * Rv); - if (Details::PhaseUsed::oil(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s) const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs); @@ -772,14 +433,14 @@ namespace Opm { voidage_rates[io] = surface_rates[io]; - if (Details::PhaseUsed::gas(pu)) { + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { voidage_rates[io] -= Rv * surface_rates[ig]; } voidage_rates[io] /= den; } - if (Details::PhaseUsed::gas(pu)) { + if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s) const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv); @@ -787,7 +448,7 @@ namespace Opm { voidage_rates[ig] = surface_rates[ig]; - if (Details::PhaseUsed::oil(pu)) { + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { voidage_rates[ig] -= Rs * surface_rates[io]; } @@ -852,7 +513,7 @@ namespace Opm { double saltConcentration; }; - Details::RegionAttributes attr_; + RegionAttributeHelpers::RegionAttributes attr_; }; } // namespace RateConverter diff --git a/opm/simulators/wells/RegionAttributeHelpers.hpp b/opm/simulators/wells/RegionAttributeHelpers.hpp new file mode 100644 index 000000000..b7d2d67a5 --- /dev/null +++ b/opm/simulators/wells/RegionAttributeHelpers.hpp @@ -0,0 +1,410 @@ +/* + Copyright 2014, 2015 SINTEF ICT, Applied Mathematics. + Copyright 2014, 2015 Statoil ASA. + Copyright 2017, IRIS + Copyright 2017, Equinor + + 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_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED +#define OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm { + namespace RegionAttributeHelpers { + /** + * Convenience tools for processing region + * spesific attributes + */ + namespace Select { + template + struct RegionIDParameter + { + using type = + typename std::remove_reference::type &; + }; + + template + struct RegionIDParameter + { + using type = RegionID; + }; + } // Select + + /** + * \brief Computes the temperature, pressure, and counter increment. + * + * In a parallel run only cells owned contribute to the cell average. + * \tparam is_parallel Whether this is a parallel run. + */ + template + struct AverageIncrementCalculator + { + /** + * \brief Computes the temperature, pressure, and counter increment. + * \param pressure The pressure. + * \param temperature The temperature. + * \param rs The rs. + * \param rv The rv. + * \param cell The current cell index. + * \param ownership A vector indicating whether a cell is owned + * by this process (value 1), or not (value 0). + * \param cell The cell index. + */ + std::tuple + operator()(const std::vector& pressure, + const std::vector& temperature, + const std::vector& rs, + const std::vector& rv, + const std::vector& ownership, + std::size_t cell){ + if ( ownership[cell] ) + { + return std::make_tuple(pressure[cell], + temperature[cell], + rs[cell], + rv[cell], + 1); + } + else + { + return std::make_tuple(0, 0, 0, 0, 0); + } + } + }; + template<> + struct AverageIncrementCalculator + { + std::tuple + operator()(const std::vector& pressure, + const std::vector& temperature, + const std::vector& rs, + const std::vector& rv, + const std::vector&, + std::size_t cell){ + return std::make_tuple(pressure[cell], + temperature[cell], + rs[cell], + rv[cell], + 1); + } + }; + /** + * Provide mapping from Region IDs to user-specified collection + * of per-region attributes. + * + * \tparam RegionId Region identifier type. Must be hashable by + * \code std::hash<> \endcode. Typically a built-in integer + * type--e.g., \c int. + * + * \tparam Attributes User-defined type that represents + * collection of attributes that have meaning in a per-region + * aggregate sense. Must be copy-constructible. + */ + template + class RegionAttributes + { + public: + /** + * Expose \c RegionId as a vocabulary type for use in query + * methods. + */ + using RegionID = + typename Select::RegionIDParameter + ::value>::type; + + using ID = + typename std::remove_reference::type; + + /** + * Aggregate per-region attributes along with region's + * representative cell. + */ + struct Value { + Value(const Attributes& attr) + : attr_(attr) + , cell_(-1) + {} + + Attributes attr_; + int cell_; + }; + + using AttributeMap = + std::unordered_map>; + + + /** + * Constructor. + * + * \tparam RMap Class type that implements the RegionMapping + * protocol. Typically an instantiation of \code + * Opm::RegionMapping<> \endcode. + * + * \param[in] rmap Specific region mapping that provides + * reverse lookup from regions to cells. + * + * \param[in] attr Pre-constructed initialiser for \c + * Attributes. + */ + template + RegionAttributes(const RMap& rmap, + const Attributes& attr) + { + using VT = typename AttributeMap::value_type; + + for (const auto& r : rmap.activeRegions()) { + auto v = std::make_unique(attr); + + const auto stat = attr_.insert(VT(r, std::move(v))); + + if (stat.second) { + // New value inserted. + const auto& cells = rmap.cells(r); + + assert (! cells.empty()); + + // Region's representative cell. + stat.first->second->cell_ = cells[0]; + } + } + } + + /** + * Retrieve representative cell in region. + * + * \param[in] reg Specific region. + * + * \return Representative cell in region \p reg. + */ + int cell(const RegionID reg) const + { + return this->find(reg).cell_; + } + + bool has(const RegionID reg) const + { + return this->attr_.find(reg) != this->attr_.end(); + } + + void insert(const RegionID r, const Attributes& attr) + { + auto [pos, inserted] = this->attr_.try_emplace(r, std::make_unique(attr)); + if (inserted) { + pos->second->cell_ = -1; // NOT -1.0 -- "cell_" is 'int' + } + } + + /** + * Request read-only access to region's attributes. + * + * + * \return Read-only access to all regions attributes. + */ + const AttributeMap& attributes() const + { + return attr_; + } + + + /** + * Request read-only access to region's attributes. + * + * \param[in] reg Specific region. + * + * \return Read-only access to region \p reg's per-region + * attributes. + */ + const Attributes& attributes(const RegionID reg) const + { + return this->find(reg).attr_; + } + + /** + * Request modifiable access to region's attributes. + * + * \param[in] reg Specific region. + * + * \return Read-write access to region \p reg's per-region + * attributes. + */ + Attributes& attributes(const RegionID reg) + { + return this->find(reg).attr_; + } + + private: + + AttributeMap attr_; + + /** + * Read-only access to region's properties. + */ + const Value& find(const RegionID reg) const + { + const auto& i = attr_.find(reg); + + if (i == attr_.end()) { + throw std::invalid_argument("Unknown region ID"); + } + + return *i->second; + } + + /** + * Read-write access to region's properties. + */ + Value& find(const RegionID reg) + { + const auto& i = attr_.find(reg); + + if (i == attr_.end()) { + throw std::invalid_argument("Unknown region ID"); + } + + return *i->second; + } + }; + + /** + * Convenience functions for querying presence/absence of + * active phases. + */ + namespace PhaseUsed { + /** + * Active water predicate. + * + * \param[in] pu Active phase object. + * + * \return Whether or not water is an active phase. + */ + inline bool + water(const PhaseUsage& pu) + { + return pu.phase_used[ BlackoilPhases::Aqua ] != 0; + } + + /** + * Active oil predicate. + * + * \param[in] pu Active phase object. + * + * \return Whether or not oil is an active phase. + */ + inline bool + oil(const PhaseUsage& pu) + { + return pu.phase_used[ BlackoilPhases::Liquid ] != 0; + } + + /** + * Active gas predicate. + * + * \param[in] pu Active phase object. + * + * \return Whether or not gas is an active phase. + */ + inline bool + gas(const PhaseUsage& pu) + { + return pu.phase_used[ BlackoilPhases::Vapour ] != 0; + } + } // namespace PhaseUsed + + /** + * Convenience functions for querying numerical IDs + * ("positions") of active phases. + */ + namespace PhasePos { + /** + * Numerical ID of active water phase. + * + * \param[in] pu Active phase object. + * + * \return Non-negative index/position of water if + * active, -1 if not. + */ + inline int + water(const PhaseUsage& pu) + { + int p = -1; + + if (PhaseUsed::water(pu)) { + p = pu.phase_pos[ BlackoilPhases::Aqua ]; + } + + return p; + } + + /** + * Numerical ID of active oil phase. + * + * \param[in] pu Active phase object. + * + * \return Non-negative index/position of oil if + * active, -1 if not. + */ + inline int + oil(const PhaseUsage& pu) + { + int p = -1; + + if (PhaseUsed::oil(pu)) { + p = pu.phase_pos[ BlackoilPhases::Liquid ]; + } + + return p; + } + + /** + * Numerical ID of active gas phase. + * + * \param[in] pu Active phase object. + * + * \return Non-negative index/position of gas if + * active, -1 if not. + */ + inline int + gas(const PhaseUsage& pu) + { + int p = -1; + + if (PhaseUsed::gas(pu)) { + p = pu.phase_pos[ BlackoilPhases::Vapour ]; + } + + return p; + } + } // namespace PhasePos + + } // namespace RegionAttributesHelpers +} // namespace Opm + +#endif /* OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED */ diff --git a/opm/simulators/wells/RegionAverageCalculator.hpp b/opm/simulators/wells/RegionAverageCalculator.hpp new file mode 100644 index 000000000..430811bbb --- /dev/null +++ b/opm/simulators/wells/RegionAverageCalculator.hpp @@ -0,0 +1,260 @@ +/* + Copyright 2021, Equinor + + 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_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED +#define OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +/** + * \file + * Facility for converting component rates at surface conditions to + * phase (voidage) rates at reservoir conditions. + * + * This uses the average hydrocarbon pressure to define fluid + * properties. The facility is intended to support Reservoir Voidage + * rates only ('RESV'). + */ + +namespace Opm { + namespace RegionAverageCalculator { + + /** + * Computes hydrocarbon weighed average pressures over regions + * + * \tparam FluidSystem Fluid system class. Expected to be a BlackOilFluidSystem + * + * \tparam Region Type of a forward region mapping. Expected + * to provide indexed access through \code operator[]() + * \endcode as well as inner types \c value_type, \c + * size_type, and \c const_iterator. Typically \code + * std::vector \endcode. + */ + template + class AverageRegionalPressure { + public: + /** + * Constructor. + * + * \param[in] region Forward region mapping. Often + * corresponds to the "FIPNUM" mapping of an ECLIPSE input + * deck. + */ + AverageRegionalPressure(const PhaseUsage& phaseUsage, + const Region& region) + : phaseUsage_(phaseUsage) + , rmap_ (region) + , attr_ (rmap_, Attributes()) + { + } + + + /** + * Compute pore volume averaged hydrocarbon state pressure, * + */ + template + void defineState(const EbosSimulator& simulator) + { + + + int numRegions = 0; + const auto& gridView = simulator.gridView(); + const auto& comm = gridView.comm(); + for (const auto& reg : rmap_.activeRegions()) { + numRegions = std::max(numRegions, reg); + } + numRegions = comm.max(numRegions); + for (int reg = 1; reg <= numRegions ; ++ reg) { + if(!attr_.has(reg)) + attr_.insert(reg, Attributes()); + } + // create map from cell to region + // and set all attributes to zero + for (int reg = 1; reg <= numRegions ; ++ reg) { + auto& ra = attr_.attributes(reg); + ra.pressure = 0.0; + ra.pv = 0.0; + + } + + // quantities for pore volume average + std::unordered_map attributes_pv; + + // quantities for hydrocarbon volume average + std::unordered_map attributes_hpv; + + for (int reg = 1; reg <= numRegions ; ++ reg) { + attributes_pv.insert({reg, Attributes()}); + attributes_hpv.insert({reg, Attributes()}); + } + + ElementContext elemCtx( simulator ); + const auto& elemEndIt = gridView.template end(); + for (auto elemIt = gridView.template begin(); + elemIt != elemEndIt; + ++elemIt) + { + + const auto& elem = *elemIt; + if (elem.partitionType() != Dune::InteriorEntity) + continue; + + elemCtx.updatePrimaryStencil(elem); + elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0); + const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0); + const auto& fs = intQuants.fluidState(); + // use pore volume weighted averages. + const double pv_cell = + simulator.model().dofTotalVolume(cellIdx) + * intQuants.porosity().value(); + + // only count oil and gas filled parts of the domain + double hydrocarbon = 1.0; + const auto& pu = phaseUsage_; + if (RegionAttributeHelpers::PhaseUsed::water(pu)) { + hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value(); + } + + const int reg = rmap_.region(cellIdx); + assert(reg >= 0); + + // sum p, rs, rv, and T. + const double hydrocarbonPV = pv_cell*hydrocarbon; + if (hydrocarbonPV > 0.) { + auto& attr = attributes_hpv[reg]; + attr.pv += hydrocarbonPV; + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV; + } else { + assert(RegionAttributeHelpers::PhaseUsed::gas(pu)); + attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV; + } + } + + if (pv_cell > 0.) { + auto& attr = attributes_pv[reg]; + attr.pv += pv_cell; + if (RegionAttributeHelpers::PhaseUsed::oil(pu)) { + attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell; + } else if (RegionAttributeHelpers::PhaseUsed::gas(pu)) { + attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell; + } else { + assert(RegionAttributeHelpers::PhaseUsed::water(pu)); + attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell; + } + } + } + + for (int reg = 1; reg <= numRegions ; ++ reg) { + auto& ra = attr_.attributes(reg); + const double hpv_sum = comm.sum(attributes_hpv[reg].pv); + // TODO: should we have some epsilon here instead of zero? + if (hpv_sum > 0.) { + const auto& attri_hpv = attributes_hpv[reg]; + const double p_hpv_sum = comm.sum(attri_hpv.pressure); + ra.pressure = p_hpv_sum / hpv_sum; + } else { + // using the pore volume to do the averaging + const auto& attri_pv = attributes_pv[reg]; + const double pv_sum = comm.sum(attri_pv.pv); + assert(pv_sum > 0.); + const double p_pv_sum = comm.sum(attri_pv.pressure); + ra.pressure = p_pv_sum / pv_sum; + + } + } + } + + /** + * Region identifier. + * + * Integral type. + */ + typedef typename RegionMapping::RegionId RegionId; + + /** + * Average pressure + * + */ + double + pressure(const RegionId r) const + { + if (r == 0 ) // region 0 is the whole field + { + double pressure = 0.0; + int num_active_regions = 0; + for (const auto& attr : attr_.attributes()) { + const auto& value = *attr.second; + const auto& ra = value.attr_; + pressure += ra.pressure; + num_active_regions ++; + } + return pressure / num_active_regions; + } + + const auto& ra = attr_.attributes(r); + return ra.pressure; + } + + + private: + /** + * Fluid property object. + */ + const PhaseUsage phaseUsage_; + + /** + * "Fluid-in-place" region mapping (forward and reverse). + */ + const RegionMapping rmap_; + + /** + * Derived property attributes for each active region. + */ + struct Attributes { + Attributes() + : pressure(0.0) + , pv(0.0) + + {} + + double pressure; + double pv; + + }; + + RegionAttributeHelpers::RegionAttributes attr_; + + }; + } // namespace RegionAverageCalculator +} // namespace Opm + +#endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */ diff --git a/opm/simulators/wells/TargetCalculator.cpp b/opm/simulators/wells/TargetCalculator.cpp index 7a5904cf0..ac52d88c5 100644 --- a/opm/simulators/wells/TargetCalculator.cpp +++ b/opm/simulators/wells/TargetCalculator.cpp @@ -40,13 +40,19 @@ namespace WellGroupHelpers TargetCalculator::TargetCalculator(const Group::ProductionCMode cmode, const PhaseUsage& pu, const std::vector& resv_coeff, - const double group_grat_target_from_sales) + const double group_grat_target_from_sales, + const std::string& group_name, + const GroupState& group_state, + const bool use_gpmaint) : cmode_(cmode) , pu_(pu) , resv_coeff_(resv_coeff) , group_grat_target_from_sales_(group_grat_target_from_sales) -{ + , group_name_(group_name) + , group_state_(group_state) + , use_gpmaint_(use_gpmaint) +{ } template @@ -107,7 +113,12 @@ double TargetCalculator::groupTarget(const Group::ProductionControls ctrl) const case Group::ProductionCMode::LRAT: return ctrl.liquid_target; case Group::ProductionCMode::RESV: + { + if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_)) + return this->group_state_.gpmaint_target(this->group_name_); + return ctrl.resv_target; + } default: // Should never be here. assert(false); @@ -142,6 +153,7 @@ InjectionTargetCalculator::InjectionTargetCalculator(const Group::InjectionCMode const double sales_target, const GroupState& group_state, const Phase& injection_phase, + const bool use_gpmaint, DeferredLogger& deferred_logger) : cmode_(cmode) , pu_(pu) @@ -149,6 +161,8 @@ InjectionTargetCalculator::InjectionTargetCalculator(const Group::InjectionCMode , group_name_(group_name) , sales_target_(sales_target) , group_state_(group_state) + , use_gpmaint_(use_gpmaint) + { // initialize to avoid warning pos_ = pu.phase_pos[BlackoilPhases::Aqua]; @@ -182,8 +196,14 @@ double InjectionTargetCalculator::groupTarget(const Group::InjectionControls& ct { switch (cmode_) { case Group::InjectionCMode::RATE: + if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_)) + return this->group_state_.gpmaint_target(this->group_name_); + return ctrl.surface_max_rate; case Group::InjectionCMode::RESV: + if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_)) + return this->group_state_.gpmaint_target(this->group_name_) / resv_coeff_[pos_]; + return ctrl.resv_max_rate / resv_coeff_[pos_]; case Group::InjectionCMode::REIN: { double production_rate = this->group_state_.injection_rein_rates(ctrl.reinj_group)[pos_]; diff --git a/opm/simulators/wells/TargetCalculator.hpp b/opm/simulators/wells/TargetCalculator.hpp index fe2566759..40b8bbd3c 100644 --- a/opm/simulators/wells/TargetCalculator.hpp +++ b/opm/simulators/wells/TargetCalculator.hpp @@ -45,7 +45,10 @@ namespace WellGroupHelpers TargetCalculator(const Group::ProductionCMode cmode, const PhaseUsage& pu, const std::vector& resv_coeff, - const double group_grat_target_from_sales); + const double group_grat_target_from_sales, + const std::string& group_name, + const GroupState& group_state, + const bool use_gpmaint); template RateType calcModeRateFromRates(const std::vector& rates) const @@ -65,6 +68,9 @@ namespace WellGroupHelpers const PhaseUsage& pu_; const std::vector& resv_coeff_; const double group_grat_target_from_sales_; + const std::string& group_name_; + const GroupState& group_state_; + bool use_gpmaint_; }; /// Based on a group control mode, extract or calculate rates, and @@ -79,6 +85,7 @@ namespace WellGroupHelpers const double sales_target, const GroupState& group_state, const Phase& injection_phase, + const bool use_gpmaint, DeferredLogger& deferred_logger); template @@ -98,6 +105,7 @@ namespace WellGroupHelpers const std::string& group_name_; double sales_target_; const GroupState& group_state_; + bool use_gpmaint_; int pos_; GuideRateModel::Target target_; }; diff --git a/opm/simulators/wells/WellGroupHelpers.cpp b/opm/simulators/wells/WellGroupHelpers.cpp index b40c7957c..c935dee14 100644 --- a/opm/simulators/wells/WellGroupHelpers.cpp +++ b/opm/simulators/wells/WellGroupHelpers.cpp @@ -21,13 +21,10 @@ #include #include -#include #include #include -#include #include #include -#include #include #include #include @@ -163,6 +160,17 @@ namespace WellGroupHelpers 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); } @@ -546,6 +554,30 @@ namespace WellGroupHelpers group_state.update_injection_reservoir_rates(group.name(), resv); } + void updateSurfaceRatesInjectionGroups(const Group& group, + const Schedule& schedule, + const int reportStepIdx, + const WellState& wellStateNupcol, + GroupState& group_state) + { + for (const std::string& groupName : group.groups()) { + const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx); + updateSurfaceRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellStateNupcol, group_state); + } + const int np = wellStateNupcol.numPhases(); + std::vector rates(np, 0.0); + for (int phase = 0; phase < np; ++phase) { + rates[phase] = sumWellPhaseRates(false, + group, + schedule, + wellStateNupcol, + reportStepIdx, + phase, + /*isInjector*/ true); + } + group_state.update_injection_surface_rates(group.name(), rates); + } + void updateWellRates(const Group& group, const Schedule& schedule, const int reportStepIdx, @@ -1118,7 +1150,7 @@ namespace WellGroupHelpers if (group_state.has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); - TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales); + TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl)); FractionCalculator fcalc(schedule, wellState, group_state, reportStepIdx, guideRate, tcalc.guideTargetMode(), pu, true, Phase::OIL); auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); }; @@ -1245,7 +1277,7 @@ namespace WellGroupHelpers const auto& gconsale = schedule[reportStepIdx].gconsale().get(group.name(), summaryState); sales_target = gconsale.sales_target; } - InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger); + InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger); FractionCalculator fcalc(schedule, wellState, group_state, reportStepIdx, guideRate, tcalc.guideTargetMode(), pu, false, injectionPhase); auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); }; diff --git a/opm/simulators/wells/WellGroupHelpers.hpp b/opm/simulators/wells/WellGroupHelpers.hpp index 2aea55564..82e60d2f7 100644 --- a/opm/simulators/wells/WellGroupHelpers.hpp +++ b/opm/simulators/wells/WellGroupHelpers.hpp @@ -22,6 +22,11 @@ #define OPM_WELLGROUPHELPERS_HEADER_INCLUDED #include +#include +#include +#include +#include +#include #include #include @@ -151,6 +156,12 @@ namespace WellGroupHelpers WellState& wellState, GroupState& group_state); + void updateSurfaceRatesInjectionGroups(const Group& group, + const Schedule& schedule, + const int reportStepIdx, + const WellState& wellStateNupcol, + GroupState& group_state); + void updateWellRates(const Group& group, const Schedule& schedule, const int reportStepIdx, @@ -181,6 +192,84 @@ namespace WellGroupHelpers WellState& wellState, GroupState& group_state); + template + void updateGpMaintTargetForGroups(const Group& group, + const Schedule& schedule, + const RegionalValues& regional_values, + const int reportStepIdx, + const double dt, + const WellState& well_state, + GroupState& group_state) + { + 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 [name, number] = *gpm->region(); + const double error = gpm->pressure_target() - regional_values.pressure(number); + double current_rate = 0.0; + const auto& pu = well_state.phaseUsage(); + double sign = 1.0; + switch (gpm->flow_target()) { + case Opm::GPMaint::FlowTarget::RESV_PROD: + { + current_rate = -group_state.injection_vrep_rate(group.name()); + sign = -1.0; + break; + } + case Opm::GPMaint::FlowTarget::RESV_OINJ: + { + if (pu.phase_used[BlackoilPhases::Liquid]) + current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Liquid]]; + + break; + } + case Opm::GPMaint::FlowTarget::RESV_WINJ: + { + if (pu.phase_used[BlackoilPhases::Aqua]) + current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Aqua]]; + + break; + } + case Opm::GPMaint::FlowTarget::RESV_GINJ: + { + if (pu.phase_used[BlackoilPhases::Vapour]) + current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Vapour]]; + break; + } + case Opm::GPMaint::FlowTarget::SURF_OINJ: + { + if (pu.phase_used[BlackoilPhases::Liquid]) + current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Liquid]]; + + break; + } + case Opm::GPMaint::FlowTarget::SURF_WINJ: + { + if (pu.phase_used[BlackoilPhases::Aqua]) + current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Aqua]]; + + break; + } + case Opm::GPMaint::FlowTarget::SURF_GINJ: + { + if (pu.phase_used[BlackoilPhases::Vapour]) + current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Vapour]]; + + break; + } + default: + throw std::invalid_argument("Invalid Flow target type in GPMAINT"); + } + + auto& gpmaint_state = group_state.gpmaint(group.name()); + double rate = gpm->rate(gpmaint_state, current_rate, error, dt); + group_state.update_gpmaint_target(group.name(), std::max(0.0, sign * rate)); + } std::map computeNetworkPressures(const Opm::Network::ExtNetwork& network, const WellState& well_state, diff --git a/opm/simulators/wells/WellInterfaceEval.cpp b/opm/simulators/wells/WellInterfaceEval.cpp index 72b89e2f7..037d83862 100644 --- a/opm/simulators/wells/WellInterfaceEval.cpp +++ b/opm/simulators/wells/WellInterfaceEval.cpp @@ -136,7 +136,7 @@ getGroupInjectionControl(const Group& group, const auto& gconsale = schedule[baseif_.currentStep()].gconsale().get(group.name(), summaryState); sales_target = gconsale.sales_target; } - WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger); + WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, false, injectionPhase); auto localFraction = [&](const std::string& child) { @@ -230,7 +230,7 @@ getGroupProductionControl(const Group& group, if (group_state.has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); - WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales); + WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl)); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, true, Phase::OIL); auto localFraction = [&](const std::string& child) { diff --git a/opm/simulators/wells/WellInterfaceFluidSystem.cpp b/opm/simulators/wells/WellInterfaceFluidSystem.cpp index 7e70b3bee..ce7e24d2d 100644 --- a/opm/simulators/wells/WellInterfaceFluidSystem.cpp +++ b/opm/simulators/wells/WellInterfaceFluidSystem.cpp @@ -1002,7 +1002,7 @@ getGroupInjectionTargetRate(const Group& group, const auto& gconsale = schedule[currentStep()].gconsale().get(group.name(), summaryState); sales_target = gconsale.sales_target; } - WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger); + WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, currentStep(), guideRate(), tcalc.guideTargetMode(), pu, false, injectionPhase); auto localFraction = [&](const std::string& child) { @@ -1072,7 +1072,7 @@ getGroupProductionTargetRate(const Group& group, if (group_state.has_grat_sales_target(group.name())) gratTargetFromSales = group_state.grat_sales_target(group.name()); - WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales); + WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl)); WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, currentStep(), guideRate(), tcalc.guideTargetMode(), pu, true, Phase::OIL); auto localFraction = [&](const std::string& child) {