Add support for gpmaint

2026-08-26 10:47:19 -05:00 · 2021-09-07 08:52:52 +02:00
parent 87b618e93a
commit 9d2f26f7e8
12 changed files with 487 additions and 19 deletions
@@ -132,6 +132,10 @@ namespace Opm {
            using RateConverterType = RateConverter::
                SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;

+            // For computing average pressured used by gpmaint
+            using AverageRegionalPressureType = RateConverter::
+                AverageRegionalPressure<FluidSystem, std::vector<int> >;
+
            BlackoilWellModel(Simulator& ebosSimulator);

            void init();
@@ -317,6 +321,8 @@ namespace Opm {
            bool alternative_well_rate_init_{};

            std::unique_ptr<RateConverterType> rateConverter_{};
+            std::unique_ptr<AverageRegionalPressureType> gpmaint_{};
+

            SimulatorReportSingle last_report_{};

@@ -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, well_state, 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
@@ -207,6 +207,13 @@ namespace Opm {
                                                    std::vector<int>(local_num_cells_, 0)));
        rateConverter_->template defineState<ElementContext>(ebosSimulator_);

+        // Compute reservoir volumes for RESV controls.
+        const auto& fp = this->eclState_.fieldProps();
+        const auto fipnum = fp.get_int("FIPNUM");
+        gpmaint_.reset(new AverageRegionalPressureType (phase_usage_,
+                                                            fipnum));
+        gpmaint_->template defineState<ElementContext>(ebosSimulator_);
+
        {
            const auto& sched_state = this->schedule()[timeStepIdx];
            // update VFP properties
@@ -458,6 +465,12 @@ namespace Opm {
        // update the rate converter with current averages pressures etc in
        rateConverter_->template defineState<ElementContext>(ebosSimulator_);

+        gpmaint_->template defineState<ElementContext>(ebosSimulator_);
+        const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx);
+        WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup,
+                                                       schedule_, *gpmaint_, reportStepIdx, dt, this->wellState(), this->groupState());
+
+
        // calculate the well potentials
        try {
            updateWellPotentials(reportStepIdx,
@@ -472,7 +485,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);

@@ -107,6 +107,27 @@ const std::vector<double>& 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<double>& rates) {
+    if (rates.size() != this->num_phases)
+        throw std::logic_error("Wrong number of phases");
+
+    this->inj_surface_rates[gname] = rates;
+}
+
+const std::vector<double>& 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 {
@@ -52,6 +52,11 @@ public:
    void update_injection_reservoir_rates(const std::string& gname, const std::vector<double>& rates);
    const std::vector<double>& 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<double>& rates);
+    const std::vector<double>& injection_surface_rates(const std::string& gname) const;
+
+
    void update_injection_rein_rates(const std::string& gname, const std::vector<double>& rates);
    const std::vector<double>& 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<std::string, Group::ProductionCMode> production_controls;
    std::map<std::string, std::vector<double>> prod_red_rates;
    std::map<std::string, std::vector<double>> inj_red_rates;
+    std::map<std::string, std::vector<double>> inj_surface_rates;
    std::map<std::string, std::vector<double>> inj_resv_rates;
    std::map<std::string, std::vector<double>> inj_potentials;
    std::map<std::string, std::vector<double>> inj_rein_rates;
    std::map<std::string, double> inj_vrep_rate;
    std::map<std::string, double> m_grat_sales_target;
+    std::map<std::string, double> m_gpmaint_target;


    std::map<std::pair<Phase, std::string>, Group::InjectionCMode> injection_controls;
@@ -197,6 +197,32 @@ namespace Opm {
                    return this->find(reg).cell_;
                }

+                bool has(const RegionID reg) const
+                {
+                    const auto& i = attr_.find(reg);
+
+                    if (i == attr_.end()) {
+                        return false;
+                    }
+
+                    return true;
+                }
+
+                void insert(const RegionID r, const Attributes& attr)
+                {
+                    using VT = typename AttributeMap::value_type;
+
+                    auto v = std::make_unique<Value>(attr);
+
+                    const auto stat = attr_.insert(VT(r, std::move(v)));
+
+                    if (stat.second) {
+                        // Region's representative cell.
+                        stat.first->second->cell_ = -1.0;
+                    }
+                }
+
+
                /**
                 * Request read-only access to region's attributes.
                 *
@@ -855,6 +881,213 @@ namespace Opm {
            Details::RegionAttributes<RegionId, Attributes> attr_;

        };
+
+        /**
+         * 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<int> \endcode.
+         */
+        template <class FluidSystem, class Region>
+        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 <typename ElementContext, class EbosSimulator>
+            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 = 0; 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 = 0; reg < numRegions; ++ reg) {
+                    auto& ra = attr_.attributes(reg);
+                    ra.pressure = 0.0;
+                    ra.pv = 0.0;
+
+                }
+
+                // quantities for pore volume average
+                std::unordered_map<RegionId, Attributes> attributes_pv;
+
+                // quantities for hydrocarbon volume average
+                std::unordered_map<RegionId, Attributes> attributes_hpv;
+
+                for (int reg = 0; reg < numRegions; ++ reg) {
+                    attributes_pv.insert({reg, Attributes()});
+                    attributes_hpv.insert({reg, Attributes()});
+                }
+
+                for (int reg = 0; reg < numRegions; ++ reg) {
+                    auto& ra = attributes_pv[reg];
+                    ra.pressure = 0.0;
+                    ra.pv = 0.0;
+                }
+                for (int reg = 0; reg < numRegions; ++ reg) {
+                    auto& ra = attributes_hpv[reg];
+                    ra.pressure = 0.0;
+                    ra.pv = 0.0;
+                }
+
+                ElementContext elemCtx( simulator );
+                 const auto& elemEndIt = gridView.template end</*codim=*/0>();
+                for (auto elemIt = gridView.template begin</*codim=*/0>();
+                     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 (Details::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 (Details::PhaseUsed::oil(pu)) {
+                            attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
+                        } else {
+                            assert(Details::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 (Details::PhaseUsed::oil(pu)) {
+                            attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
+                        } else if (Details::PhaseUsed::gas(pu)) {
+                             attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
+                        } else {
+                            assert(Details::PhaseUsed::water(pu));
+                            attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
+                        }
+                    }
+                }
+
+                for (int reg = 0; 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<Region>::RegionId RegionId;
+
+            /**
+             * Average pressure
+             *
+             */
+            double
+            fpr(const RegionId r) const
+            {
+                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<Region> rmap_;
+
+            /**
+             * Derived property attributes for each active region.
+             */
+            struct Attributes {
+                Attributes()
+                    : pressure   (0.0)
+                    , pv(0.0)
+
+                {}
+
+                double pressure;
+                double pv;
+
+            };
+
+            Details::RegionAttributes<RegionId, Attributes> attr_;
+
+        };
    } // namespace RateConverter
 } // namespace Opm

@@ -40,13 +40,19 @@ namespace WellGroupHelpers
 TargetCalculator::TargetCalculator(const Group::ProductionCMode cmode,
                                   const PhaseUsage& pu,
                                   const std::vector<double>& 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 <typename RateType>
@@ -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_];
@@ -45,7 +45,10 @@ namespace WellGroupHelpers
        TargetCalculator(const Group::ProductionCMode cmode,
                         const PhaseUsage& pu,
                         const std::vector<double>& 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 <typename RateType>
        RateType calcModeRateFromRates(const std::vector<RateType>& rates) const
@@ -65,6 +68,9 @@ namespace WellGroupHelpers
        const PhaseUsage& pu_;
        const std::vector<double>& 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 <typename RateVec>
@@ -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_;
    };
@@ -21,13 +21,10 @@
 #include <config.h>
 #include <opm/simulators/wells/WellGroupHelpers.hpp>

-#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Group/GConSump.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Group/GConSale.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/Group/Group.hpp>
 #include <opm/simulators/utils/DeferredLogger.hpp>
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
-#include <opm/simulators/wells/GroupState.hpp>
 #include <opm/simulators/wells/TargetCalculator.hpp>
 #include <opm/simulators/wells/VFPProdProperties.hpp>
 #include <opm/simulators/wells/WellState.hpp>
@@ -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,31 @@ 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,
+                                             WellState& wellState,
+                                             GroupState& group_state)
+    {
+        for (const std::string& groupName : group.groups()) {
+            const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
+            updateSurfaceRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellStateNupcol, wellState, group_state);
+        }
+        const int np = wellState.numPhases();
+        std::vector<double> 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 +1151,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 +1278,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); };
@@ -22,6 +22,11 @@
 #define OPM_WELLGROUPHELPERS_HEADER_INCLUDED

 #include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group/GPMaint.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/simulators/wells/GroupState.hpp>
+#include <opm/simulators/wells/WellState.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Group/Group.hpp>

 #include <map>
 #include <string>
@@ -151,6 +156,13 @@ namespace WellGroupHelpers
                                             WellState& wellState,
                                             GroupState& group_state);

+    void updateSurfaceRatesInjectionGroups(const Group& group,
+                                             const Schedule& schedule,
+                                             const int reportStepIdx,
+                                             const WellState& wellStateNupcol,
+                                             WellState& wellState,
+                                             GroupState& group_state);
+
    void updateWellRates(const Group& group,
                         const Schedule& schedule,
                         const int reportStepIdx,
@@ -181,6 +193,84 @@ namespace WellGroupHelpers
                             WellState& wellState,
                             GroupState& group_state);

+    template <class RegionalValues>
+    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.fpr(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<std::string, double>
    computeNetworkPressures(const Opm::Network::ExtNetwork& network,
                            const WellState& well_state,
@@ -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) {
@@ -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) {