add WellInterfaceEval

CMakeLists_files.cmake

@@ -66,6 +66,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/VFPProdProperties.cpp
   opm/simulators/wells/VFPInjProperties.cpp
   opm/simulators/wells/WellGroupHelpers.cpp
+  opm/simulators/wells/WellInterfaceEval.cpp
   opm/simulators/wells/WellInterfaceFluidSystem.cpp
   opm/simulators/wells/WellInterfaceGeneric.cpp
   opm/simulators/wells/WellInterfaceIndices.cpp

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -25,6 +25,8 @@
 
 #include <opm/parser/eclipse/Units/UnitSystem.hpp>
 
+#include <opm/simulators/wells/VFPProperties.hpp>
+
 #include <algorithm>
 #include <utility>
 

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -1789,7 +1789,7 @@ namespace Opm
             // Setup function for evaluation of BHP from THP (used only if needed).
             auto bhp_from_thp = [&]() {
                 const auto rates = getRates();
-                return calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
+                return this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
             };
             // Call generic implementation.
             Base::assembleControlEqInj(well_state, group_state, schedule, summaryState, inj_controls, getBhp(), injection_rate, bhp_from_thp, control_eq, deferred_logger);
@@ -1798,7 +1798,7 @@ namespace Opm
             const auto rates = getRates();
             // Setup function for evaluation of BHP from THP (used only if needed).
             auto bhp_from_thp = [&]() {
-                return calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
+                return this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
             };
             // Call generic implementation.
             Base::assembleControlEqProd(well_state, group_state, schedule, summaryState, prod_controls, getBhp(), rates, bhp_from_thp, control_eq, deferred_logger);

opm/simulators/wells/StandardWell.hpp

@@ -31,6 +31,7 @@
 #include <opm/simulators/wells/RateConverter.hpp>
 #include <opm/simulators/wells/VFPInjProperties.hpp>
 #include <opm/simulators/wells/VFPProdProperties.hpp>
+#include <opm/simulators/wells/WellHelpers.hpp>
 #include <opm/simulators/wells/WellInterface.hpp>
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>

opm/simulators/wells/StandardWell_impl.hpp

@@ -893,7 +893,7 @@ namespace Opm
             // Setup function for evaluation of BHP from THP (used only if needed).
             auto bhp_from_thp = [&]() {
                 const auto rates = getRates();
-                return calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
+                return this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
             };
             // Call generic implementation.
             const auto& inj_controls = well.injectionControls(summaryState);
@@ -903,7 +903,7 @@ namespace Opm
             const auto rates = getRates();
             // Setup function for evaluation of BHP from THP (used only if needed).
             auto bhp_from_thp = [&]() {
-                 return calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
+                 return this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
             };
             // Call generic implementation.
             const auto& prod_controls = well.productionControls(summaryState);

opm/simulators/wells/WellInterface.hpp

@@ -30,13 +30,9 @@
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTestState.hpp>
-#include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp>
 
 #include <opm/core/props/BlackoilPhases.hpp>
 
-#include <opm/simulators/wells/VFPProperties.hpp>
-#include <opm/simulators/wells/WellHelpers.hpp>
-#include <opm/simulators/wells/WellGroupHelpers.hpp>
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/WellState.hpp>
 // NOTE: GasLiftSingleWell.hpp includes StandardWell.hpp which includes ourself
@@ -56,15 +52,11 @@ namespace Opm {
 #include<dune/istl/bcrsmatrix.hh>
 #include<dune/istl/matrixmatrix.hh>
 
-#include <opm/material/densead/Math.hpp>
 #include <opm/material/densead/Evaluation.hpp>
 
 #include <opm/simulators/wells/WellInterfaceIndices.hpp>
 
-#include <array>
 #include <cassert>
-#include <memory>
-#include <string>
 #include <vector>
 
 namespace Opm
@@ -232,7 +224,6 @@ public:
 
     Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;
 
-
     template <class EvalWell>
     Eval restrictEval(const EvalWell& in) const
     {
@@ -300,9 +291,6 @@ protected:
 
     double wsalt() const;
 
-    template <class ValueType>
-    ValueType calculateBhpFromThp(const WellState& well_state, const std::vector<ValueType>& rates, const Well& well, const SummaryState& summaryState, DeferredLogger& deferred_logger) const;
-
     virtual double getRefDensity() const = 0;
 
     // Component fractions for each phase for the well
@@ -355,54 +343,6 @@ protected:
 
     void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
                              DeferredLogger& deferred_logger);
-
-    template <class EvalWell>
-    void getGroupInjectionControl(const Group& group,
-                                  const WellState& well_state,
-                                  const GroupState& group_state,
-                                  const Schedule& schedule,
-                                  const SummaryState& summaryState,
-                                  const InjectorType& injectorType,
-                                  const EvalWell& bhp,
-                                  const EvalWell& injection_rate,
-                                  EvalWell& control_eq,
-                                  double efficiencyFactor,
-                                  DeferredLogger& deferred_logger);
-
-    template <class EvalWell>
-    void getGroupProductionControl(const Group& group,
-                                   const WellState& well_state,
-                                   const GroupState& group_state,
-                                   const Schedule& schedule,
-                                   const SummaryState& summaryState,
-                                   const EvalWell& bhp,
-                                   const std::vector<EvalWell>& rates,
-                                   EvalWell& control_eq,
-                                   double efficiencyFactor);
-
-    template <class EvalWell, class BhpFromThpFunc>
-    void assembleControlEqInj(const WellState& well_state,
-                              const GroupState& group_state,
-                              const Schedule& schedule,
-                              const SummaryState& summaryState,
-                              const Well::InjectionControls& controls,
-                              const EvalWell& bhp,
-                              const EvalWell& injection_rate,
-                              BhpFromThpFunc bhp_from_thp,
-                              EvalWell& control_eq,
-                              DeferredLogger& deferred_logger);
-
-    template <class EvalWell, class BhpFromThpFunc>
-    void assembleControlEqProd(const WellState& well_state,
-                               const GroupState& group_state,
-                               const Schedule& schedule,
-                               const SummaryState& summaryState,
-                               const Well::ProductionControls& controls,
-                               const EvalWell& bhp,
-                               const std::vector<EvalWell>& rates, // Always 3 canonical rates.
-                               BhpFromThpFunc bhp_from_thp,
-                               EvalWell& control_eq,
-                               DeferredLogger& deferred_logger);
 };
 
 }

opm/simulators/wells/WellInterfaceEval.cpp

@@ -0,0 +1,572 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2018 IRIS
+
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+#include <opm/simulators/wells/WellInterfaceEval.hpp>
+
+#include <ebos/eclalternativeblackoilindices.hh>
+
+#include <opm/material/densead/Evaluation.hpp>
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+
+#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
+#include <opm/simulators/wells/GroupState.hpp>
+#include <opm/simulators/wells/RateConverter.hpp>
+#include <opm/simulators/wells/TargetCalculator.hpp>
+#include <opm/simulators/wells/VFPProperties.hpp>
+#include <opm/simulators/wells/WellGroupHelpers.hpp>
+#include <opm/simulators/wells/WellHelpers.hpp>
+#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
+#include <opm/simulators/wells/WellState.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <stdexcept>
+
+namespace Opm
+{
+
+template<class FluidSystem>
+WellInterfaceEval<FluidSystem>::
+WellInterfaceEval(const WellInterfaceFluidSystem<FluidSystem>& baseif)
+    : baseif_(baseif)
+{}
+
+template<class FluidSystem>
+template<class EvalWell>
+void
+WellInterfaceEval<FluidSystem>::
+getGroupInjectionControl(const Group& group,
+                         const WellState& well_state,
+                         const GroupState& group_state,
+                         const Schedule& schedule,
+                         const SummaryState& summaryState,
+                         const InjectorType& injectorType,
+                         const EvalWell& bhp,
+                         const EvalWell& injection_rate,
+                         EvalWell& control_eq,
+                         double efficiencyFactor,
+                         DeferredLogger& deferred_logger)
+{
+    // Setting some defaults to silence warnings below.
+    // Will be overwritten in the switch statement.
+    Phase injectionPhase = Phase::WATER;
+    switch (injectorType) {
+    case InjectorType::WATER:
+    {
+        injectionPhase = Phase::WATER;
+        break;
+    }
+    case InjectorType::OIL:
+    {
+        injectionPhase = Phase::OIL;
+        break;
+    }
+    case InjectorType::GAS:
+    {
+        injectionPhase = Phase::GAS;
+        break;
+    }
+    default:
+        // Should not be here.
+        assert(false);
+    }
+
+    auto currentGroupControl = group_state.injection_control(group.name(), injectionPhase);
+    if (currentGroupControl == Group::InjectionCMode::FLD ||
+        currentGroupControl == Group::InjectionCMode::NONE) {
+        if (!group.injectionGroupControlAvailable(injectionPhase)) {
+            // We cannot go any further up the hierarchy. This could
+            // be the FIELD group, or any group for which this has
+            // been set in GCONINJE or GCONPROD. If we are here
+            // anyway, it is likely that the deck set inconsistent
+            // requirements, such as GRUP control mode on a well with
+            // no appropriate controls defined on any of its
+            // containing groups. We will therefore use the wells' bhp
+            // limit equation as a fallback.
+            const auto& controls = baseif_.wellEcl().injectionControls(summaryState);
+            control_eq = bhp - controls.bhp_limit;
+            return;
+        } else {
+            // Inject share of parents control
+            const auto& parent = schedule.getGroup( group.parent(), baseif_.currentStep());
+            efficiencyFactor *= group.getGroupEfficiencyFactor();
+            getGroupInjectionControl(parent, well_state, group_state, schedule, summaryState, injectorType, bhp, injection_rate, control_eq, efficiencyFactor, deferred_logger);
+            return;
+        }
+    }
+
+    efficiencyFactor *= group.getGroupEfficiencyFactor();
+    const auto& well = baseif_.wellEcl();
+    const auto pu = baseif_.phaseUsage();
+
+    if (!group.isInjectionGroup()) {
+        // use bhp as control eq and let the updateControl code find a valid control
+        const auto& controls = well.injectionControls(summaryState);
+        control_eq = bhp - controls.bhp_limit;
+        return;
+    }
+
+    // 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.
+
+    // Make conversion factors for RESV <-> surface rates.
+    std::vector<double> resv_coeff(pu.num_phases, 1.0);
+    baseif_.rateConverter().calcCoeff(0, baseif_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+
+    double sales_target = 0;
+    if (schedule[baseif_.currentStep()].gconsale().has(group.name())) {
+        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::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, false, injectionPhase);
+
+    auto localFraction = [&](const std::string& child) {
+        return fcalc.localFraction(child, "");
+    };
+
+    auto localReduction = [&](const std::string& group_name) {
+        const std::vector<double>& groupTargetReductions = group_state.injection_reduction_rates(group_name);
+        return tcalc.calcModeRateFromRates(groupTargetReductions);
+    };
+
+    const double orig_target = tcalc.groupTarget(group.injectionControls(injectionPhase, summaryState), deferred_logger);
+    const auto chain = WellGroupHelpers::groupChainTopBot(baseif_.name(), group.name(), schedule, baseif_.currentStep());
+    // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
+    const size_t num_ancestors = chain.size() - 1;
+    double target = orig_target;
+    for (size_t ii = 0; ii < num_ancestors; ++ii) {
+        if ((ii == 0) || baseif_.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.
+            target -= localReduction(chain[ii]);
+        }
+        target *= localFraction(chain[ii+1]);
+    }
+    // Avoid negative target rates coming from too large local reductions.
+    const double target_rate = std::max(0.0, target / efficiencyFactor);
+    const auto current_rate = injection_rate; // Switch sign since 'rates' are negative for producers.
+    control_eq = current_rate - target_rate;
+}
+
+template<class FluidSystem>
+template<class EvalWell>
+void
+WellInterfaceEval<FluidSystem>::
+getGroupProductionControl(const Group& group,
+                          const WellState& well_state,
+                          const GroupState& group_state,
+                          const Schedule& schedule,
+                          const SummaryState& summaryState,
+                          const EvalWell& bhp,
+                          const std::vector<EvalWell>& rates,
+                          EvalWell& control_eq,
+                          double efficiencyFactor)
+{
+    const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name());
+    if (currentGroupControl == Group::ProductionCMode::FLD ||
+        currentGroupControl == Group::ProductionCMode::NONE) {
+        if (!group.productionGroupControlAvailable()) {
+            // We cannot go any further up the hierarchy. This could
+            // be the FIELD group, or any group for which this has
+            // been set in GCONINJE or GCONPROD. If we are here
+            // anyway, it is likely that the deck set inconsistent
+            // requirements, such as GRUP control mode on a well with
+            // no appropriate controls defined on any of its
+            // containing groups. We will therefore use the wells' bhp
+            // limit equation as a fallback.
+            const auto& controls = baseif_.wellEcl().productionControls(summaryState);
+            control_eq = bhp - controls.bhp_limit;
+            return;
+        } else {
+            // Produce share of parents control
+            const auto& parent = schedule.getGroup(group.parent(), baseif_.currentStep());
+            efficiencyFactor *= group.getGroupEfficiencyFactor();
+            getGroupProductionControl(parent, well_state, group_state, schedule, summaryState, bhp, rates, control_eq, efficiencyFactor);
+            return;
+        }
+    }
+
+    efficiencyFactor *= group.getGroupEfficiencyFactor();
+    const auto& well = baseif_.wellEcl();
+    const auto pu = baseif_.phaseUsage();
+
+    if (!group.isProductionGroup()) {
+        // use bhp as control eq and let the updateControl code find a valid control
+        const auto& controls = well.productionControls(summaryState);
+        control_eq = bhp - controls.bhp_limit;
+        return;
+    }
+
+    // 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.
+
+    // Make conversion factors for RESV <-> surface rates.
+    std::vector<double> resv_coeff(baseif_.phaseUsage().num_phases, 1.0);
+    baseif_.rateConverter().calcCoeff(0, baseif_.pvtRegionIdx(), resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
+
+    // gconsale may adjust the grat target.
+    // the adjusted rates is send to the targetCalculator
+    double gratTargetFromSales = 0.0;
+    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::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, true, Phase::OIL);
+
+    auto localFraction = [&](const std::string& child) {
+        return fcalc.localFraction(child, "");
+    };
+
+    auto localReduction = [&](const std::string& group_name) {
+        const std::vector<double>& groupTargetReductions = group_state.production_reduction_rates(group_name);
+        return tcalc.calcModeRateFromRates(groupTargetReductions);
+    };
+
+    const double orig_target = tcalc.groupTarget(group.productionControls(summaryState));
+    const auto chain = WellGroupHelpers::groupChainTopBot(baseif_.name(), group.name(), schedule, baseif_.currentStep());
+    // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
+    const size_t num_ancestors = chain.size() - 1;
+    double target = orig_target;
+    for (size_t ii = 0; ii < num_ancestors; ++ii) {
+        if ((ii == 0) || baseif_.guideRate()->has(chain[ii])) {
+            // Apply local reductions only at the control level
+            // (top) and for levels where we have a specified
+            // group guide rate.
+            target -= localReduction(chain[ii]);
+        }
+        target *= localFraction(chain[ii+1]);
+    }
+    // Avoid negative target rates coming from too large local reductions.
+    const double target_rate = std::max(0.0, target / efficiencyFactor);
+    const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
+    control_eq = current_rate - target_rate;
+}
+
+template<class FluidSystem>
+template<class EvalWell>
+void
+WellInterfaceEval<FluidSystem>::
+assembleControlEqProd_(const WellState& well_state,
+                       const GroupState& group_state,
+                       const Schedule& schedule,
+                       const SummaryState& summaryState,
+                       const Well::ProductionControls& controls,
+                       const EvalWell& bhp,
+                       const std::vector<EvalWell>& rates, // Always 3 canonical rates.
+                       const std::function<EvalWell()>& bhp_from_thp,
+                       EvalWell& control_eq,
+                       DeferredLogger& deferred_logger)
+{
+    auto current = well_state.currentProductionControl(baseif_.indexOfWell());
+    const auto& pu = baseif_.phaseUsage();
+    const double efficiencyFactor = baseif_.wellEcl().getEfficiencyFactor();
+
+    switch (current) {
+    case Well::ProducerCMode::ORAT: {
+        assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
+        const EvalWell rate = -rates[BlackoilPhases::Liquid];
+        control_eq = rate - controls.oil_rate;
+        break;
+    }
+    case Well::ProducerCMode::WRAT: {
+        assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
+        const EvalWell rate = -rates[BlackoilPhases::Aqua];
+        control_eq = rate - controls.water_rate;
+        break;
+    }
+    case Well::ProducerCMode::GRAT: {
+        assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
+        const EvalWell rate = -rates[BlackoilPhases::Vapour];
+        control_eq = rate - controls.gas_rate;
+        break;
+    }
+    case Well::ProducerCMode::LRAT: {
+        assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
+        assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
+        EvalWell rate = -rates[BlackoilPhases::Aqua] - rates[BlackoilPhases::Liquid];
+        control_eq = rate - controls.liquid_rate;
+        break;
+    }
+    case Well::ProducerCMode::CRAT: {
+        OPM_DEFLOG_THROW(std::runtime_error, "CRAT control not supported " << baseif_.name(), deferred_logger);
+    }
+    case Well::ProducerCMode::RESV: {
+        auto total_rate = rates[0]; // To get the correct type only.
+        total_rate = 0.0;
+        std::vector<double> convert_coeff(baseif_.numPhases(), 1.0);
+        baseif_.rateConverter().calcCoeff(/*fipreg*/ 0, baseif_.pvtRegionIdx(), convert_coeff);
+        for (int phase = 0; phase < 3; ++phase) {
+            if (pu.phase_used[phase]) {
+                const int pos = pu.phase_pos[phase];
+                total_rate -= rates[phase] * convert_coeff[pos]; // Note different indices.
+            }
+        }
+        if (controls.prediction_mode) {
+            control_eq = total_rate - controls.resv_rate;
+        } else {
+            std::vector<double> hrates(baseif_.numPhases(), 0.);
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+                hrates[pu.phase_pos[Water]] = controls.water_rate;
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                hrates[pu.phase_pos[Oil]] = controls.oil_rate;
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                hrates[pu.phase_pos[Gas]] = controls.gas_rate;
+            }
+            std::vector<double> hrates_resv(baseif_.numPhases(), 0.);
+            baseif_.rateConverter().calcReservoirVoidageRates(/*fipreg*/ 0, baseif_.pvtRegionIdx(), hrates, hrates_resv);
+            double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
+            control_eq = total_rate - target;
+        }
+        break;
+    }
+    case Well::ProducerCMode::BHP: {
+        control_eq = bhp - controls.bhp_limit;
+        break;
+    }
+    case Well::ProducerCMode::THP: {
+        control_eq = bhp - bhp_from_thp();
+        break;
+    }
+    case Well::ProducerCMode::GRUP: {
+        assert(baseif_.wellEcl().isAvailableForGroupControl());
+        const auto& group = schedule.getGroup(baseif_.wellEcl().groupName(), baseif_.currentStep());
+        // Annoying thing: the rates passed to this function are
+        // always of size 3 and in canonical (for PhaseUsage)
+        // order. This is what is needed for VFP calculations if
+        // they are required (THP controlled well). But for the
+        // group production control things we must pass only the
+        // active phases' rates.
+        std::vector<EvalWell> active_rates(pu.num_phases);
+        for (int canonical_phase = 0; canonical_phase < 3; ++canonical_phase) {
+            if (pu.phase_used[canonical_phase]) {
+                active_rates[pu.phase_pos[canonical_phase]] = rates[canonical_phase];
+            }
+        }
+        getGroupProductionControl(group, well_state, group_state, schedule, summaryState, bhp, active_rates, control_eq, efficiencyFactor);
+        break;
+    }
+    case Well::ProducerCMode::CMODE_UNDEFINED: {
+        OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + baseif_.name(), deferred_logger);
+    }
+    case Well::ProducerCMode::NONE: {
+        OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + baseif_.name(), deferred_logger);
+    }
+    }
+}
+
+template<class FluidSystem>
+template<class EvalWell>
+void
+WellInterfaceEval<FluidSystem>::
+assembleControlEqInj_(const WellState& well_state,
+                      const GroupState& group_state,
+                      const Schedule& schedule,
+                      const SummaryState& summaryState,
+                      const Well::InjectionControls& controls,
+                      const EvalWell& bhp,
+                      const EvalWell& injection_rate,
+                      const std::function<EvalWell()>& bhp_from_thp,
+                      EvalWell& control_eq,
+                      DeferredLogger& deferred_logger)
+{
+    auto current = well_state.currentInjectionControl(baseif_.indexOfWell());
+    const InjectorType injectorType = controls.injector_type;
+    const auto& pu = baseif_.phaseUsage();
+    const double efficiencyFactor = baseif_.wellEcl().getEfficiencyFactor();
+
+    switch (current) {
+    case Well::InjectorCMode::RATE: {
+        control_eq = injection_rate - controls.surface_rate;
+        break;
+    }
+    case Well::InjectorCMode::RESV: {
+        std::vector<double> convert_coeff(baseif_.numPhases(), 1.0);
+        baseif_.rateConverter().calcCoeff(/*fipreg*/ 0, baseif_.pvtRegionIdx(), convert_coeff);
+
+        double coeff = 1.0;
+
+        switch (injectorType) {
+        case InjectorType::WATER: {
+            coeff = convert_coeff[pu.phase_pos[BlackoilPhases::Aqua]];
+            break;
+        }
+        case InjectorType::OIL: {
+            coeff = convert_coeff[pu.phase_pos[BlackoilPhases::Liquid]];
+            break;
+        }
+        case InjectorType::GAS: {
+            coeff = convert_coeff[pu.phase_pos[BlackoilPhases::Vapour]];
+            break;
+        }
+        default:
+            throw("Expected WATER, OIL or GAS as type for injectors " + baseif_.wellEcl().name());
+        }
+
+        control_eq = coeff * injection_rate - controls.reservoir_rate;
+        break;
+    }
+    case Well::InjectorCMode::THP: {
+        control_eq = bhp - bhp_from_thp();
+        break;
+    }
+    case Well::InjectorCMode::BHP: {
+        control_eq = bhp - controls.bhp_limit;
+        break;
+    }
+    case Well::InjectorCMode::GRUP: {
+        assert(baseif_.wellEcl().isAvailableForGroupControl());
+        const auto& group = schedule.getGroup(baseif_.wellEcl().groupName(), baseif_.currentStep());
+        this->getGroupInjectionControl(group,
+                                       well_state,
+                                       group_state,
+                                       schedule,
+                                       summaryState,
+                                       injectorType,
+                                       bhp,
+                                       injection_rate,
+                                       control_eq,
+                                       efficiencyFactor,
+                                       deferred_logger);
+        break;
+    }
+    case Well::InjectorCMode::CMODE_UNDEFINED: {
+        OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + baseif_.name(), deferred_logger);
+    }
+    }
+}
+
+template<class FluidSystem>
+template<class EvalWell>
+EvalWell
+WellInterfaceEval<FluidSystem>::
+calculateBhpFromThp(const WellState& well_state,
+                    const std::vector<EvalWell>& rates,
+                    const Well& well,
+                    const SummaryState& summaryState,
+                    const double rho,
+                    DeferredLogger& deferred_logger) const
+{
+    // TODO: when well is under THP control, the BHP is dependent on the rates,
+    // the well rates is also dependent on the BHP, so it might need to do some iteration.
+    // However, when group control is involved, change of the rates might impacts other wells
+    // so iterations on a higher level will be required. Some investigation might be needed when
+    // we face problems under THP control.
+
+    assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
+
+    const EvalWell aqua = rates[Water];
+    const EvalWell liquid = rates[Oil];
+    const EvalWell vapour = rates[Gas];
+
+    // pick the reference density
+    // typically the reference in the top layer
+    if (baseif_.isInjector() )
+    {
+        const auto& controls = well.injectionControls(summaryState);
+        const double vfp_ref_depth = baseif_.vfpProperties()->getInj()->getTable(controls.vfp_table_number).getDatumDepth();
+        const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
+        return baseif_.vfpProperties()->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, baseif_.getTHPConstraint(summaryState)) - dp;
+     }
+     else if (baseif_.isProducer()) {
+         const auto& controls = well.productionControls(summaryState);
+         const double vfp_ref_depth = baseif_.vfpProperties()->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
+         const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
+         return baseif_.vfpProperties()->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, baseif_.getTHPConstraint(summaryState), baseif_.getALQ(well_state)) - dp;
+     }
+     else {
+         OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + baseif_.name(), deferred_logger);
+     }
+}
+
+#define INSTANCE_METHODS(A,...) \
+template void WellInterfaceEval<A>:: \
+assembleControlEqProd_<__VA_ARGS__>(const WellState&, \
+                                    const GroupState&, \
+                                    const Schedule&, \
+                                    const SummaryState&, \
+                                    const Well::ProductionControls&, \
+                                    const __VA_ARGS__&, \
+                                    const std::vector<__VA_ARGS__>&, \
+                                    const std::function<__VA_ARGS__()>&, \
+                                    __VA_ARGS__&, \
+                                    DeferredLogger&); \
+template void WellInterfaceEval<A>:: \
+assembleControlEqInj_<__VA_ARGS__>(const WellState&, \
+                                   const GroupState&, \
+                                   const Schedule&, \
+                                   const SummaryState&, \
+                                   const Well::InjectionControls&, \
+                                   const __VA_ARGS__&, \
+                                   const __VA_ARGS__&, \
+                                   const std::function<__VA_ARGS__()>&, \
+                                   __VA_ARGS__&, \
+                                   DeferredLogger&); \
+template __VA_ARGS__ WellInterfaceEval<A>:: \
+calculateBhpFromThp<__VA_ARGS__>(const WellState&, \
+                                 const std::vector<__VA_ARGS__>&, \
+                                 const Well&, \
+                                 const SummaryState&, \
+                                 const double, \
+                                 DeferredLogger&) const;
+
+using FluidSys = BlackOilFluidSystem<double, BlackOilDefaultIndexTraits>;
+using FluidAltSys = BlackOilFluidSystem<double, EclAlternativeBlackOilIndexTraits>;
+
+template class WellInterfaceEval<FluidSys>;
+template class WellInterfaceEval<FluidAltSys>;
+
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,3,0u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,4,0u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,5,0u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,6,0u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,7,0u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,8,0u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,4u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,5u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,6u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,7u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,8u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,9u>)
+INSTANCE_METHODS(FluidSys, DenseAd::Evaluation<double,-1,10u>)
+
+INSTANCE_METHODS(FluidAltSys, DenseAd::Evaluation<double,7,0u>)
+INSTANCE_METHODS(FluidAltSys, DenseAd::Evaluation<double,-1,8u>)
+
+#define INSTANCE_BHP(...) \
+template double WellInterfaceEval<__VA_ARGS__>:: \
+calculateBhpFromThp<double>(const WellState&, \
+                            const std::vector<double>&, \
+                            const Well&, \
+                            const SummaryState&, \
+                            const double, \
+                            DeferredLogger&) const;
+
+INSTANCE_BHP(FluidSys)
+INSTANCE_BHP(FluidAltSys)
+
+} // namespace Opm

opm/simulators/wells/WellInterfaceEval.hpp

@@ -0,0 +1,166 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2017 IRIS
+  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 <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPM_WELLINTERFACE_EVAL_HEADER_INCLUDED
+#define OPM_WELLINTERFACE_EVAL_HEADER_INCLUDED
+
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleTypes.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
+
+#include <functional>
+
+namespace Opm
+{
+
+class DeferredLogger;
+class Group;
+class GroupState;
+class Schedule;
+class SummaryState;
+template<class FluidSystem> class WellInterfaceFluidSystem;
+class WellState;
+
+template<class FluidSystem>
+class WellInterfaceEval {
+    static constexpr int Water = BlackoilPhases::Aqua;
+    static constexpr int Oil = BlackoilPhases::Liquid;
+    static constexpr int Gas = BlackoilPhases::Vapour;
+
+protected:
+    WellInterfaceEval(const WellInterfaceFluidSystem<FluidSystem>& baseif);
+
+    template<class EvalWell>
+    void getGroupInjectionControl(const Group& group,
+                                  const WellState& well_state,
+                                  const GroupState& group_state,
+                                  const Schedule& schedule,
+                                  const SummaryState& summaryState,
+                                  const InjectorType& injectorType,
+                                  const EvalWell& bhp,
+                                  const EvalWell& injection_rate,
+                                  EvalWell& control_eq,
+                                  double efficiencyFactor,
+                                  DeferredLogger& deferred_logger);
+
+
+    template<class EvalWell>
+    void getGroupProductionControl(const Group& group,
+                                   const WellState& well_state,
+                                   const GroupState& group_state,
+                                   const Schedule& schedule,
+                                   const SummaryState& summaryState,
+                                   const EvalWell& bhp,
+                                   const std::vector<EvalWell>& rates,
+                                   EvalWell& control_eq,
+                                   double efficiencyFactor);
+
+    template<class EvalWell, class BhpFromThpFunc>
+    void assembleControlEqProd(const WellState& well_state,
+                               const GroupState& group_state,
+                               const Schedule& schedule,
+                               const SummaryState& summaryState,
+                               const Well::ProductionControls& controls,
+                               const EvalWell& bhp,
+                               const std::vector<EvalWell>& rates, // Always 3 canonical rates.
+                               BhpFromThpFunc bhp_from_thp,
+                               EvalWell& control_eq,
+                               DeferredLogger& deferred_logger)
+    {
+        std::function<EvalWell()> eval = [&bhp_from_thp]() { return bhp_from_thp(); };
+        assembleControlEqProd_(well_state,
+                               group_state,
+                               schedule,
+                               summaryState,
+                               controls,
+                               bhp,
+                               rates,
+                               eval,
+                               control_eq,
+                               deferred_logger);
+    }
+
+    template<class EvalWell>
+    void assembleControlEqProd_(const WellState& well_state,
+                                const GroupState& group_state,
+                                const Schedule& schedule,
+                                const SummaryState& summaryState,
+                                const Well::ProductionControls& controls,
+                                const EvalWell& bhp,
+                                const std::vector<EvalWell>& rates, // Always 3 canonical rates.
+                                const std::function<EvalWell()>& bhp_from_thp,
+                                EvalWell& control_eq,
+                                DeferredLogger& deferred_logger);
+
+    template<class EvalWell, class BhpFromThpFunc>
+    void assembleControlEqInj(const WellState& well_state,
+                              const GroupState& group_state,
+                              const Schedule& schedule,
+                              const SummaryState& summaryState,
+                              const Well::InjectionControls& controls,
+                              const EvalWell& bhp,
+                              const EvalWell& injection_rate,
+                              BhpFromThpFunc bhp_from_thp,
+                              EvalWell& control_eq,
+                              DeferredLogger& deferred_logger)
+    {
+        std::function<EvalWell()> eval = [&bhp_from_thp]() { return bhp_from_thp(); };
+        assembleControlEqInj_(well_state,
+                              group_state,
+                              schedule,
+                              summaryState,
+                              controls,
+                              bhp,
+                              injection_rate,
+                              eval,
+                              control_eq,
+                              deferred_logger);
+    }
+
+    template<class EvalWell>
+    void assembleControlEqInj_(const WellState& well_state,
+                               const GroupState& group_state,
+                               const Schedule& schedule,
+                               const SummaryState& summaryState,
+                               const Well::InjectionControls& controls,
+                               const EvalWell& bhp,
+                               const EvalWell& injection_rate,
+                               const std::function<EvalWell()>& bhp_from_thp,
+                               EvalWell& control_eq,
+                               DeferredLogger& deferred_logger);
+
+    template <class EvalWell>
+    EvalWell calculateBhpFromThp(const WellState& well_state,
+                                 const std::vector<EvalWell>& rates,
+                                 const Well& well,
+                                 const SummaryState& summaryState,
+                                 const double rho,
+                                 DeferredLogger& deferred_logger) const;
+
+    const WellInterfaceFluidSystem<FluidSystem>& baseif_;
+};
+
+}
+
+#endif // OPM_WELLINTERFACE_EVAL_HEADER_INCLUDED

opm/simulators/wells/WellInterfaceFluidSystem.hpp

@@ -43,6 +43,12 @@ class WellState;
 
 template<class FluidSystem>
 class WellInterfaceFluidSystem : public WellInterfaceGeneric {
+protected:
+    using RateConverterType = RateConverter::
+    SurfaceToReservoirVoidage<FluidSystem, std::vector<int>>;
+    // to indicate a invalid completion
+    static constexpr int INVALIDCOMPLETION = std::numeric_limits<int>::max();
+
 public:
     void updateWellTestState(const WellState& well_state,
                              const double& simulationTime,
@@ -56,13 +62,12 @@ public:
     static constexpr int Oil = BlackoilPhases::Liquid;
     static constexpr int Gas = BlackoilPhases::Vapour;
 
+    const RateConverterType& rateConverter() const
+    {
+        return rateConverter_;
+    }
+
 protected:
-    using RateConverterType = RateConverter::
-    SurfaceToReservoirVoidage<FluidSystem, std::vector<int>>;
-
-    // to indicate a invalid completion
-    static constexpr int INVALIDCOMPLETION = std::numeric_limits<int>::max();
-
     WellInterfaceFluidSystem(const Well& well,
                              const ParallelWellInfo& parallel_well_info,
                              const int time_step,

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -107,6 +107,37 @@ public:
         return this->wellStatus_ == Well::Status::STOP;
     }
 
+    int currentStep() const {
+        return this->current_step_;
+    }
+
+    int pvtRegionIdx() const {
+        return pvtRegionIdx_;
+    }
+
+    const GuideRate* guideRate() const {
+        return guide_rate_;
+    }
+
+    int numPhases() const {
+        return number_of_phases_;
+    }
+
+    double refDepth() const {
+        return ref_depth_;
+    }
+
+    const VFPProperties* vfpProperties() const {
+        return vfp_properties_;
+    }
+
+    double gravity() const {
+        return gravity_;
+    }
+
+    double getTHPConstraint(const SummaryState& summaryState) const;
+    double getALQ(const WellState& well_state) const;
+
 protected:
     // whether a well is specified with a non-zero and valid VFP table number
     bool isVFPActive(DeferredLogger& deferred_logger) const;
@@ -114,13 +145,11 @@ protected:
     bool getAllowCrossFlow() const;
     double wsolvent() const;
     double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
-    double getTHPConstraint(const SummaryState& summaryState) const;
     void updateWellTestStatePhysical(const WellState& well_state,
                                      const double simulation_time,
                                      const bool write_message_to_opmlog,
                                      WellTestState& well_test_state,
                                      DeferredLogger& deferred_logger) const;
-    double getALQ(const WellState& well_state) const;
 
     // definition of the struct OperabilityStatus
     struct OperabilityStatus {

opm/simulators/wells/WellInterfaceIndices.cpp

@@ -25,6 +25,7 @@
 #include <ebos/eclalternativeblackoilindices.hh>
 
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+
 #include <opm/models/blackoil/blackoilindices.hh>
 #include <opm/models/blackoil/blackoilonephaseindices.hh>
 #include <opm/models/blackoil/blackoiltwophaseindices.hh>
@@ -56,6 +57,7 @@ WellInterfaceIndices(const Well& well,
                                             index_of_well,
                                             first_perf_index,
                                             perf_data)
+    , WellInterfaceEval<FluidSystem>(static_cast<const WellInterfaceFluidSystem<FluidSystem>&>(*this))
 {
 }
 

opm/simulators/wells/WellInterfaceIndices.hpp

@@ -24,12 +24,15 @@
 #define OPM_WELLINTERFACE_INDICES_HEADER_INCLUDED
 
 #include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
+#include <opm/simulators/wells/WellInterfaceEval.hpp>
 
 namespace Opm
 {
 
 template<class FluidSystem, class Indices, class Scalar>
-class WellInterfaceIndices : public WellInterfaceFluidSystem<FluidSystem> {
+class WellInterfaceIndices : public WellInterfaceFluidSystem<FluidSystem>
+                           , public WellInterfaceEval<FluidSystem>
+{
 public:
     using WellInterfaceFluidSystem<FluidSystem>::Gas;
     using WellInterfaceFluidSystem<FluidSystem>::Oil;
@@ -51,9 +54,6 @@ protected:
     int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
 
     double scalingFactor(const int phaseIdx) const;
-
-    Scalar volumetricSurfaceRateForConnection(int cellIdx,
-                                              int phaseIdx) const;
 };
 
 }

opm/simulators/wells/WellInterface_impl.hpp

@@ -211,54 +211,6 @@ namespace Opm
 
 
 
-    template<typename TypeTag>
-    template<class ValueType>
-    ValueType
-    WellInterface<TypeTag>::
-    calculateBhpFromThp(const WellState &well_state,
-                        const std::vector<ValueType>& rates,
-                        const Well& well,
-                        const SummaryState& summaryState,
-                        DeferredLogger& deferred_logger) const
-    {
-        // TODO: when well is under THP control, the BHP is dependent on the rates,
-        // the well rates is also dependent on the BHP, so it might need to do some iteration.
-        // However, when group control is involved, change of the rates might impacts other wells
-        // so iterations on a higher level will be required. Some investigation might be needed when
-        // we face problems under THP control.
-
-        assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
-
-        const ValueType aqua = rates[Water];
-        const ValueType liquid = rates[Oil];
-        const ValueType vapour = rates[Gas];
-
-        // pick the reference density
-        // typically the reference in the top layer
-        const double rho = getRefDensity();
-
-        if (this->isInjector() )
-        {
-            const auto& controls = well.injectionControls(summaryState);
-            const double vfp_ref_depth = this->vfp_properties_->getInj()->getTable(controls.vfp_table_number).getDatumDepth();
-            const double dp = wellhelpers::computeHydrostaticCorrection(this->ref_depth_, vfp_ref_depth, rho, this->gravity_);
-            return this->vfp_properties_->getInj()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState)) - dp;
-         }
-         else if (this->isProducer()) {
-             const auto& controls = well.productionControls(summaryState);
-             const double vfp_ref_depth = this->vfp_properties_->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
-             const double dp = wellhelpers::computeHydrostaticCorrection(this->ref_depth_, vfp_ref_depth, rho, this->gravity_);
-             return this->vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState), this->getALQ(well_state)) - dp;
-         }
-         else {
-             OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + this->name(), deferred_logger);
-         }
-
-
-
-    }
-
-
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
@@ -669,7 +621,7 @@ namespace Opm
                 for (int p = 0; p<np; ++p) {
                     rates[p] = well_state.wellRates(well_index)[p];
                 }
-                double bhp = calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
+                double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
                 well_state.update_bhp(well_index, bhp);
 
                 // if the total rates are negative or zero
@@ -879,7 +831,7 @@ namespace Opm
                 for (int p = 0; p<np; ++p) {
                     rates[p] = well_state.wellRates(well_index)[p];
                 }
-                double bhp = calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
+                double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
                 well_state.update_bhp(well_index, bhp);
 
                 // if the total rates are negative or zero
@@ -954,420 +906,6 @@ namespace Opm
 
 
 
-    template <typename TypeTag>
-    template <class EvalWell, class BhpFromThpFunc>
-    void
-    WellInterface<TypeTag>::assembleControlEqInj(const WellState& well_state,
-                                                 const GroupState& group_state,
-                                                 const Schedule& schedule,
-                                                 const SummaryState& summaryState,
-                                                 const Well::InjectionControls& controls,
-                                                 const EvalWell& bhp,
-                                                 const EvalWell& injection_rate,
-                                                 BhpFromThpFunc bhp_from_thp,
-                                                 EvalWell& control_eq,
-                                                 DeferredLogger& deferred_logger)
-    {
-        auto current = well_state.currentInjectionControl(this->index_of_well_);
-        const InjectorType injectorType = controls.injector_type;
-        const auto& pu = this->phaseUsage();
-        const double efficiencyFactor = this->well_ecl_.getEfficiencyFactor();
-
-        switch (current) {
-        case Well::InjectorCMode::RATE: {
-            control_eq = injection_rate - controls.surface_rate;
-            break;
-        }
-        case Well::InjectorCMode::RESV: {
-            std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
-            this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
-
-            double coeff = 1.0;
-
-            switch (injectorType) {
-            case InjectorType::WATER: {
-                coeff = convert_coeff[pu.phase_pos[BlackoilPhases::Aqua]];
-                break;
-            }
-            case InjectorType::OIL: {
-                coeff = convert_coeff[pu.phase_pos[BlackoilPhases::Liquid]];
-                break;
-            }
-            case InjectorType::GAS: {
-                coeff = convert_coeff[pu.phase_pos[BlackoilPhases::Vapour]];
-                break;
-            }
-            default:
-                throw("Expected WATER, OIL or GAS as type for injectors " + this->well_ecl_.name());
-            }
-
-            control_eq = coeff * injection_rate - controls.reservoir_rate;
-            break;
-        }
-        case Well::InjectorCMode::THP: {
-            control_eq = bhp - bhp_from_thp();
-            break;
-        }
-        case Well::InjectorCMode::BHP: {
-            control_eq = bhp - controls.bhp_limit;
-            break;
-        }
-        case Well::InjectorCMode::GRUP: {
-            assert(this->well_ecl_.isAvailableForGroupControl());
-            const auto& group = schedule.getGroup(this->well_ecl_.groupName(), this->current_step_);
-            getGroupInjectionControl(group,
-                                     well_state,
-                                     group_state,
-                                     schedule,
-                                     summaryState,
-                                     injectorType,
-                                     bhp,
-                                     injection_rate,
-                                     control_eq,
-                                     efficiencyFactor,
-                                     deferred_logger);
-            break;
-        }
-        case Well::InjectorCMode::CMODE_UNDEFINED: {
-            OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name(), deferred_logger);
-        }
-        }
-    }
-
-
-
-
-    template <typename TypeTag>
-    template <class EvalWell, class BhpFromThpFunc>
-    void
-    WellInterface<TypeTag>::assembleControlEqProd(const WellState& well_state,
-                                                  const GroupState& group_state,
-                                                  const Schedule& schedule,
-                                                  const SummaryState& summaryState,
-                                                  const Well::ProductionControls& controls,
-                                                  const EvalWell& bhp,
-                                                  const std::vector<EvalWell>& rates, // Always 3 canonical rates.
-                                                  BhpFromThpFunc bhp_from_thp,
-                                                  EvalWell& control_eq,
-                                                  DeferredLogger& deferred_logger)
-    {
-        auto current = well_state.currentProductionControl(this->index_of_well_);
-        const auto& pu = this->phaseUsage();
-        const double efficiencyFactor = this->well_ecl_.getEfficiencyFactor();
-
-        switch (current) {
-        case Well::ProducerCMode::ORAT: {
-            assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
-            const EvalWell rate = -rates[BlackoilPhases::Liquid];
-            control_eq = rate - controls.oil_rate;
-            break;
-        }
-        case Well::ProducerCMode::WRAT: {
-            assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
-            const EvalWell rate = -rates[BlackoilPhases::Aqua];
-            control_eq = rate - controls.water_rate;
-            break;
-        }
-        case Well::ProducerCMode::GRAT: {
-            assert(FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx));
-            const EvalWell rate = -rates[BlackoilPhases::Vapour];
-            control_eq = rate - controls.gas_rate;
-            break;
-        }
-        case Well::ProducerCMode::LRAT: {
-            assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx));
-            assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx));
-            EvalWell rate = -rates[BlackoilPhases::Aqua] - rates[BlackoilPhases::Liquid];
-            control_eq = rate - controls.liquid_rate;
-            break;
-        }
-        case Well::ProducerCMode::CRAT: {
-            OPM_DEFLOG_THROW(std::runtime_error, "CRAT control not supported " << this->name(), deferred_logger);
-        }
-        case Well::ProducerCMode::RESV: {
-            auto total_rate = rates[0]; // To get the correct type only.
-            total_rate = 0.0;
-            std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
-            this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
-            for (int phase = 0; phase < 3; ++phase) {
-                if (pu.phase_used[phase]) {
-                    const int pos = pu.phase_pos[phase];
-                    total_rate -= rates[phase] * convert_coeff[pos]; // Note different indices.
-                }
-            }
-            if (controls.prediction_mode) {
-                control_eq = total_rate - controls.resv_rate;
-            } else {
-                std::vector<double> hrates(this->number_of_phases_, 0.);
-                if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                    hrates[pu.phase_pos[Water]] = controls.water_rate;
-                }
-                if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
-                    hrates[pu.phase_pos[Oil]] = controls.oil_rate;
-                }
-                if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                    hrates[pu.phase_pos[Gas]] = controls.gas_rate;
-                }
-                std::vector<double> hrates_resv(this->number_of_phases_, 0.);
-                this->rateConverter_.calcReservoirVoidageRates(/*fipreg*/ 0, this->pvtRegionIdx_, hrates, hrates_resv);
-                double target = std::accumulate(hrates_resv.begin(), hrates_resv.end(), 0.0);
-                control_eq = total_rate - target;
-            }
-            break;
-        }
-        case Well::ProducerCMode::BHP: {
-            control_eq = bhp - controls.bhp_limit;
-            break;
-        }
-        case Well::ProducerCMode::THP: {
-            control_eq = bhp - bhp_from_thp();
-            break;
-        }
-        case Well::ProducerCMode::GRUP: {
-            assert(this->well_ecl_.isAvailableForGroupControl());
-            const auto& group = schedule.getGroup(this->well_ecl_.groupName(), this->current_step_);
-            // Annoying thing: the rates passed to this function are
-            // always of size 3 and in canonical (for PhaseUsage)
-            // order. This is what is needed for VFP calculations if
-            // they are required (THP controlled well). But for the
-            // group production control things we must pass only the
-            // active phases' rates.
-            std::vector<EvalWell> active_rates(pu.num_phases);
-            for (int canonical_phase = 0; canonical_phase < 3; ++canonical_phase) {
-                if (pu.phase_used[canonical_phase]) {
-                    active_rates[pu.phase_pos[canonical_phase]] = rates[canonical_phase];
-                }
-            }
-            getGroupProductionControl(group, well_state, group_state, schedule, summaryState, bhp, active_rates, control_eq, efficiencyFactor);
-            break;
-        }
-        case Well::ProducerCMode::CMODE_UNDEFINED: {
-            OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name(), deferred_logger);
-        }
-        case Well::ProducerCMode::NONE: {
-            OPM_DEFLOG_THROW(std::runtime_error, "Well control must be specified for well " + this->name(), deferred_logger);
-        }
-        }
-    }
-
-
-
-    template <typename TypeTag>
-    template <class EvalWell>
-    void
-    WellInterface<TypeTag>::getGroupInjectionControl(const Group& group,
-                                                      const WellState& well_state,
-                                                     const GroupState& group_state,
-                                                      const Schedule& schedule,
-                                                      const SummaryState& summaryState,
-                                                      const InjectorType& injectorType,
-                                                      const EvalWell& bhp,
-                                                      const EvalWell& injection_rate,
-                                                      EvalWell& control_eq,
-                                                      double efficiencyFactor,
-                                                      DeferredLogger& deferred_logger)
-    {
-        // Setting some defaults to silence warnings below.
-        // Will be overwritten in the switch statement.
-        Phase injectionPhase = Phase::WATER;
-        switch (injectorType) {
-        case InjectorType::WATER:
-        {
-            injectionPhase = Phase::WATER;
-            break;
-        }
-        case InjectorType::OIL:
-        {
-            injectionPhase = Phase::OIL;
-            break;
-        }
-        case InjectorType::GAS:
-        {
-            injectionPhase = Phase::GAS;
-            break;
-        }
-        default:
-            // Should not be here.
-            assert(false);
-        }
-
-        auto currentGroupControl = group_state.injection_control(group.name(), injectionPhase);
-        if (currentGroupControl == Group::InjectionCMode::FLD ||
-            currentGroupControl == Group::InjectionCMode::NONE) {
-            if (!group.injectionGroupControlAvailable(injectionPhase)) {
-                // We cannot go any further up the hierarchy. This could
-                // be the FIELD group, or any group for which this has
-                // been set in GCONINJE or GCONPROD. If we are here
-                // anyway, it is likely that the deck set inconsistent
-                // requirements, such as GRUP control mode on a well with
-                // no appropriate controls defined on any of its
-                // containing groups. We will therefore use the wells' bhp
-                // limit equation as a fallback.
-                const auto& controls = this->well_ecl_.injectionControls(summaryState);
-                control_eq = bhp - controls.bhp_limit;
-                return;
-            } else {
-                // Inject share of parents control
-                const auto& parent = schedule.getGroup( group.parent(), this->current_step_ );
-                efficiencyFactor *= group.getGroupEfficiencyFactor();
-                getGroupInjectionControl(parent, well_state, group_state, schedule, summaryState, injectorType, bhp, injection_rate, control_eq, efficiencyFactor, deferred_logger);
-                return;
-            }
-        }
-
-        efficiencyFactor *= group.getGroupEfficiencyFactor();
-        const auto& well = this->well_ecl_;
-        const auto pu = this->phaseUsage();
-
-        if (!group.isInjectionGroup()) {
-            // use bhp as control eq and let the updateControl code find a valid control
-            const auto& controls = well.injectionControls(summaryState);
-            control_eq = bhp - controls.bhp_limit;
-            return;
-        }
-
-        // 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.
-
-        // Make conversion factors for RESV <-> surface rates.
-        std::vector<double> resv_coeff(this->phaseUsage().num_phases, 1.0);
-        this->rateConverter_.calcCoeff(0, this->pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
-
-        double sales_target = 0;
-        if (schedule[this->current_step_].gconsale().has(group.name())) {
-            const auto& gconsale = schedule[this->current_step_].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::FractionCalculator fcalc(schedule, well_state, group_state, this->current_step_, this->guide_rate_, tcalc.guideTargetMode(), pu, false, injectionPhase);
-
-        auto localFraction = [&](const std::string& child) {
-            return fcalc.localFraction(child, "");
-        };
-
-        auto localReduction = [&](const std::string& group_name) {
-            const std::vector<double>& groupTargetReductions = group_state.injection_reduction_rates(group_name);
-            return tcalc.calcModeRateFromRates(groupTargetReductions);
-        };
-
-        const double orig_target = tcalc.groupTarget(group.injectionControls(injectionPhase, summaryState), deferred_logger);
-        const auto chain = WellGroupHelpers::groupChainTopBot(this->name(), group.name(), schedule, this->current_step_);
-        // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
-        const size_t num_ancestors = chain.size() - 1;
-        double target = orig_target;
-        for (size_t ii = 0; ii < num_ancestors; ++ii) {
-            if ((ii == 0) || this->guide_rate_->has(chain[ii], injectionPhase)) {
-                // Apply local reductions only at the control level
-                // (top) and for levels where we have a specified
-                // group guide rate.
-                target -= localReduction(chain[ii]);
-            }
-            target *= localFraction(chain[ii+1]);
-        }
-        // Avoid negative target rates coming from too large local reductions.
-        const double target_rate = std::max(0.0, target / efficiencyFactor);
-        const auto current_rate = injection_rate; // Switch sign since 'rates' are negative for producers.
-        control_eq = current_rate - target_rate;
-    }
-
-
-
-    template <typename TypeTag>
-    template <class EvalWell>
-    void
-    WellInterface<TypeTag>::getGroupProductionControl(const Group& group,
-                                                      const WellState& well_state,
-                                                      const GroupState& group_state,
-                                                      const Schedule& schedule,
-                                                      const SummaryState& summaryState,
-                                                      const EvalWell& bhp,
-                                                      const std::vector<EvalWell>& rates,
-                                                      EvalWell& control_eq,
-                                                      double efficiencyFactor)
-    {
-        const Group::ProductionCMode& currentGroupControl = group_state.production_control(group.name());
-        if (currentGroupControl == Group::ProductionCMode::FLD ||
-            currentGroupControl == Group::ProductionCMode::NONE) {
-            if (!group.productionGroupControlAvailable()) {
-                // We cannot go any further up the hierarchy. This could
-                // be the FIELD group, or any group for which this has
-                // been set in GCONINJE or GCONPROD. If we are here
-                // anyway, it is likely that the deck set inconsistent
-                // requirements, such as GRUP control mode on a well with
-                // no appropriate controls defined on any of its
-                // containing groups. We will therefore use the wells' bhp
-                // limit equation as a fallback.
-                const auto& controls = this->well_ecl_.productionControls(summaryState);
-                control_eq = bhp - controls.bhp_limit;
-                return;
-            } else {
-                // Produce share of parents control
-                const auto& parent = schedule.getGroup( group.parent(), this->current_step_ );
-                efficiencyFactor *= group.getGroupEfficiencyFactor();
-                getGroupProductionControl(parent, well_state, group_state, schedule, summaryState, bhp, rates, control_eq, efficiencyFactor);
-                return;
-            }
-        }
-
-        efficiencyFactor *= group.getGroupEfficiencyFactor();
-        const auto& well = this->well_ecl_;
-        const auto pu = this->phaseUsage();
-
-        if (!group.isProductionGroup()) {
-            // use bhp as control eq and let the updateControl code find a valid control
-            const auto& controls = well.productionControls(summaryState);
-            control_eq = bhp - controls.bhp_limit;
-            return;
-        }
-
-        // 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.
-
-        // Make conversion factors for RESV <-> surface rates.
-        std::vector<double> resv_coeff(this->phaseUsage().num_phases, 1.0);
-        this->rateConverter_.calcCoeff(0, this->pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
-
-        // gconsale may adjust the grat target.
-        // the adjusted rates is send to the targetCalculator
-        double gratTargetFromSales = 0.0;
-        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::FractionCalculator fcalc(schedule, well_state, group_state, this->current_step_, this->guide_rate_, tcalc.guideTargetMode(), pu, true, Phase::OIL);
-
-        auto localFraction = [&](const std::string& child) {
-            return fcalc.localFraction(child, "");
-        };
-
-        auto localReduction = [&](const std::string& group_name) {
-            const std::vector<double>& groupTargetReductions = group_state.production_reduction_rates(group_name);
-            return tcalc.calcModeRateFromRates(groupTargetReductions);
-        };
-
-        const double orig_target = tcalc.groupTarget(group.productionControls(summaryState));
-        const auto chain = WellGroupHelpers::groupChainTopBot(this->name(), group.name(), schedule, this->current_step_);
-        // Because 'name' is the last of the elements, and not an ancestor, we subtract one below.
-        const size_t num_ancestors = chain.size() - 1;
-        double target = orig_target;
-        for (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.
-                target -= localReduction(chain[ii]);
-            }
-            target *= localFraction(chain[ii+1]);
-        }
-        // Avoid negative target rates coming from too large local reductions.
-        const double target_rate = std::max(0.0, target / efficiencyFactor);
-        const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
-        control_eq = current_rate - target_rate;
-    }
-
-
-
-
     template <typename TypeTag>
     void
     WellInterface<TypeTag>::