Major restructuring to improve support for groups, multiple controls etc. Work in progress.

opm/core/WellsGroup.cpp

@@ -1,36 +1,64 @@
-/* 
+/*
- * File:   WellsGroup.cpp
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
- * Author: kjetilo
+
- * 
+  This file is part of the Open Porous Media project (OPM).
- * Created on March 27, 2012, 9:27 AM
+
- */
+  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 <opm/core/WellsGroup.hpp>
 #include <cmath>
 #include <opm/core/newwells.h>
+#include <opm/core/fluid/blackoil/phaseUsageFromDeck.hpp>
 
 namespace Opm
 {
 
+    // ==========   WellPhasesSummed methods   ===========
+
     WellPhasesSummed::WellPhasesSummed() 
-    : bhp_sum(0.0), rate_sum(0.0)
     {
-        
+        for (int i = 0; i < 3; ++i) {
+            res_inj_rates[i] = 0.0;
+            res_prod_rates[i] = 0.0;
+            surf_inj_rates[i] = 0.0;
+            surf_prod_rates[i] = 0.0;
+        }
     }
-    
+
     void WellPhasesSummed::operator+=(const WellPhasesSummed& other) 
     {
-        rate_sum += other.rate_sum;
+        for (int i = 0; i < 3; ++i) {
-        bhp_sum += other.bhp_sum;
+            res_inj_rates[i] += other.res_inj_rates[i];
+            res_prod_rates[i] += other.res_prod_rates[i];
+            surf_inj_rates[i] += other.surf_inj_rates[i];
+            surf_prod_rates[i] += other.surf_prod_rates[i];
+        }
     }
-    
+
+    // ==========   WellsGroupInterface methods   ===========
+
+
     WellsGroupInterface::WellsGroupInterface(const std::string& myname,
-            ProductionSpecification prod_spec,
+                                             const ProductionSpecification& prod_spec,
-            InjectionSpecification inje_spec)
+                                             const InjectionSpecification& inje_spec,
-    : parent_(NULL),
+                                             const PhaseUsage& phase_usage)
-    name_(myname),
+        : parent_(NULL),
-    production_specification_(prod_spec),
+          name_(myname),
-    injection_specification_(inje_spec)
+          production_specification_(prod_spec),
+          injection_specification_(inje_spec),
+          phase_usage_(phase_usage)
     {
     }
 
@@ -47,13 +75,6 @@ namespace Opm
         return name_;
     }
 
-    WellsGroup::WellsGroup(const std::string& myname,
-            ProductionSpecification prod_spec,
-            InjectionSpecification inj_spec)
-    : WellsGroupInterface(myname, prod_spec, inj_spec)
-    {
-    }
-
     bool WellsGroupInterface::isLeafNode() const
     {
         return false;
@@ -90,6 +111,62 @@ namespace Opm
         return injection_specification_;
     }
 
+    /// Calculates the correct rate for the given ProductionSpecification::ControlMode
+    double WellsGroupInterface::rateByMode(const double* res_rates, 
+                                           const double* surf_rates,
+                                           const ProductionSpecification::ControlMode mode)
+    {
+        switch (mode) {
+        case ProductionSpecification::ORAT:
+            return surf_rates[phaseUsage().phase_pos[BlackoilPhases::Liquid]];
+        case ProductionSpecification::WRAT:
+            return surf_rates[phaseUsage().phase_pos[BlackoilPhases::Aqua]];
+        case ProductionSpecification::GRAT:
+            return surf_rates[phaseUsage().phase_pos[BlackoilPhases::Vapour]];
+        case ProductionSpecification::LRAT:
+            return surf_rates[phaseUsage().phase_pos[BlackoilPhases::Liquid]] 
+                + surf_rates[phaseUsage().phase_pos[BlackoilPhases::Aqua]];
+        case ProductionSpecification::RESV:
+            {
+                double tot_rate = 0.0;
+                for (int phase = 0; phase < phaseUsage().num_phases; ++phase) {
+                    tot_rate += res_rates[phase];
+                }
+                return tot_rate;
+            }
+        default:
+            THROW("No rate associated with production control mode" << mode);
+        }
+    }
+
+    /// Calculates the correct rate for the given InjectionSpecification::ControlMode
+    double WellsGroupInterface::rateByMode(const double* res_rates, 
+                                           const double* surf_rates,
+                                           const InjectionSpecification::ControlMode mode)
+    {
+        const double* rates = 0;
+        switch (mode) {
+        case InjectionSpecification::RATE:
+            rates = surf_rates;
+            break;
+        case InjectionSpecification::RESV:
+            rates = res_rates;
+            break;
+        default:
+            THROW("No rate associated with injection control mode" << mode);
+        }
+        double tot_rate = 0.0;
+        for (int phase = 0; phase < phaseUsage().num_phases; ++phase) {
+            tot_rate += rates[phase];
+        }
+        return tot_rate;
+    }
+
+
+   
+
+    // ==============   WellsGroup members =============
+
     WellsGroupInterface* WellsGroup::findGroup(const std::string& name_of_node)
     {
         if (name() == name_of_node) {
@@ -107,89 +184,233 @@ namespace Opm
         }
     }
 
-    
+
-    void WellsGroup::calculateGuideRates()
+    WellsGroup::WellsGroup(const std::string& myname,
+                           const ProductionSpecification& prod_spec,
+                           const InjectionSpecification& inj_spec,
+                           const PhaseUsage& phase_usage)
+        : WellsGroupInterface(myname, prod_spec, inj_spec, phase_usage)
     {
-        double guide_rate_sum = 0.0;
-        for(size_t i = 0; i < children_.size(); i++) {
-            if(children_[i]->isLeafNode()) {
-                guide_rate_sum += children_[i]->prodSpec().guide_rate_;
-            }
-            else
-            {
-                children_[i]->calculateGuideRates();
-            }
-        }
-        if(guide_rate_sum != 0.0) {
-            for(size_t i = 0; i < children_.size(); i++) {
-                children_[i]->prodSpec().guide_rate_ /= guide_rate_sum;
-            }
-        }
-    }
-    
-    void WellsGroup::applyControl(const WellControlType type) 
-    {
-        for (size_t i = 0; i < children_.size(); ++i) {
-            children_[i]->applyControl(type);
-        }
     }
 
-    
+
-    bool WellsGroup::conditionsMet(const std::vector<double>& well_bhp, 
+    void WellsGroup::calculateGuideRates()
-                                   const std::vector<double>& well_rate, 
-                                   WellPhasesSummed& summed_phases, 
-                                   const double epsilon)
     {
+        double inj_guide_rate_sum = 0.0;
+        double prod_guide_rate_sum = 0.0;
+        for (size_t i = 0; i < children_.size(); i++) {
+            children_[i]->calculateGuideRates();
+            inj_guide_rate_sum += children_[i]->injSpec().guide_rate_;
+            prod_guide_rate_sum += children_[i]->prodSpec().guide_rate_;
+        }
+        injSpec().guide_rate_ = inj_guide_rate_sum;
+        prodSpec().guide_rate_ = prod_guide_rate_sum;
+    }
+    
+
+
+    /// Sets the current active control to the provided one for all injectors within the group.
+    /// After this call, the combined rate (which rate depending on control_mode) of the group
+    /// shall be equal to target.
+    void WellsGroup::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
+                                          const double target)
+    {
+        for (size_t i = 0; i < children_.size(); ++i) {
+            const double child_target = target * children_[i]->injSpec().guide_rate_/injSpec().guide_rate_;
+            children_[i]->applyInjGroupControl(control_mode, child_target);
+        }
+        injSpec().control_mode_ = InjectionSpecification::FLD;
+    }
+
+    /// Sets the current active control to the provided one for all producers within the group.
+    /// After this call, the combined rate (which rate depending on control_mode) of the group
+    /// shall be equal to target.
+    void WellsGroup::applyProdGroupControl(const ProductionSpecification::ControlMode control_mode,
+                                           const double target)
+    {
+        for (size_t i = 0; i < children_.size(); ++i) {
+            const double child_target = target * children_[i]->prodSpec().guide_rate_/prodSpec().guide_rate_;
+            children_[i]->applyProdGroupControl(control_mode, child_target);
+        }
+        prodSpec().control_mode_ = ProductionSpecification::FLD;
+    }
+
+
+    bool WellsGroup::conditionsMet(const std::vector<double>& well_bhp,
+                                   const std::vector<double>& well_reservoirrates_phase,
+                                   const std::vector<double>& well_surfacerates_phase,
+                                   WellPhasesSummed& summed_phases)
+    {
+        // Check children's constraints recursively.
         WellPhasesSummed child_phases_summed;
-        for(size_t i = 0; i < children_.size(); ++i) {
+        for (size_t i = 0; i < children_.size(); ++i) {
             WellPhasesSummed current_child_phases_summed;
-            if(!children_[i]->conditionsMet(well_bhp, well_rate, 
+            if (!children_[i]->conditionsMet(well_bhp,
-                                            current_child_phases_summed, epsilon)) {
+                                             well_reservoirrates_phase,
+                                             well_surfacerates_phase,
+                                             current_child_phases_summed)) {
                 return false;
             }
             child_phases_summed += current_child_phases_summed;
         }
-        
 
-        double bhp_target = std::min(std::abs(injSpec().BHP_limit_), prodSpec().BHP_limit_);
+        const int np = phaseUsage().num_phases;
+
+        // Injection constraints.
+        // RATE
+        if (injSpec().control_mode_ != InjectionSpecification::RATE) {
+            const double target_rate = injSpec().surface_flow_max_rate_;
+            if (target_rate >= 0.0) {
+                double my_rate = 0.0;
+                for (int phase = 0; phase < np; ++phase) {
+                    my_rate += child_phases_summed.surf_inj_rates[phase];
+                }
+                if (my_rate > target_rate) {
+                    std::cout << "Group RATE target not met for group " << name() << std::endl;
+                    std::cout << "target = " << target_rate << '\n'
+                              << "rate   = " << my_rate << std::endl;
+                    applyInjGroupControl(InjectionSpecification::RATE, target_rate);
+                    injSpec().control_mode_ = InjectionSpecification::RATE;
+                    return false;
+                }
+            }
+        }
+        // RESV
+        if (injSpec().control_mode_ != InjectionSpecification::RESV) {
+            const double target_rate = injSpec().reservoir_flow_max_rate_;
+            if (target_rate >= 0.0) {
+                double my_rate = 0.0;
+                for (int phase = 0; phase < np; ++phase) {
+                    my_rate += child_phases_summed.res_inj_rates[phase];
+                }
+                if (my_rate > target_rate) {
+                    std::cout << "Group RESV target not met for group " << name() << std::endl;
+                    std::cout << "target = " << target_rate << '\n'
+                              << "rate   = " << my_rate << std::endl;
+                    applyInjGroupControl(InjectionSpecification::RESV, target_rate);
+                    injSpec().control_mode_ = InjectionSpecification::RESV;
+                    return false;
+                }
+            }
+        }
+        // REIN
+        // \TODO: Add support for REIN controls.
+
+        // Production constraints.
+        bool prod_restrictions_violated = false;
+        ProductionSpecification::ControlMode violated_prod_mode = ProductionSpecification::NONE;
+
+        rateByMode(child_phases_summed.res_prod_rates,
+                   child_phases_summed.surf_prod_rates,
+                   mode);
+        
+        // ORAT
+        if (prodSpec().control_mode_ != ProductionSpecification::ORAT) {
+            const double target_rate = prodSpec().oil_max_rate_;
+            if (target_rate >= 0.0) {
+                const double my_rate
+                    = child_phases_summed.surf_prod_rates[phaseUsage().phase_pos[BlackoilPhases::Liquid]];
+                if (std::fabs(my_rate) > target_rate) {
+                    std::cout << "Group ORAT target not met for group " << name() << std::endl;
+                    std::cout << "target = " << target_rate << '\n'
+                              << "rate   = " << my_rate << std::endl;
+                    applyProdGroupControl(ProductionSpecification::ORAT, target_rate);
+                    prodSpec().control_mode_ = ProductionSpecification::ORAT;
+                    return false;
+                }
+            }
+        }
+
+        // WRAT
+        if (prodSpec().control_mode_ != ProductionSpecification::WRAT) {
+            const double target_rate = prodSpec().water_max_rate_;
+            if (target_rate >= 0.0) {
+                const double my_rate
+                    = child_phases_summed.surf_prod_rates[phaseUsage().phase_pos[BlackoilPhases::Aqua]];
+                if (std::fabs(my_rate) > target_rate) {
+                    std::cout << "Group WRAT target not met for group " << name() << std::endl;
+                    std::cout << "target = " << target_rate << '\n'
+                              << "rate   = " << my_rate << std::endl;
+                    applyProdGroupControl(ProductionSpecification::WRAT, target_rate);
+                    prodSpec().control_mode_ = ProductionSpecification::WRAT;
+                    return false;
+                }
+            }
+        }
+
+        // GRAT
+        if (prodSpec().control_mode_ != ProductionSpecification::GRAT) {
+            const double target_rate = prodSpec().gas_max_rate_;
+            if (target_rate >= 0.0) {
+                const double my_rate
+                    = child_phases_summed.surf_prod_rates[phaseUsage().phase_pos[BlackoilPhases::Vapour]];
+                if (std::fabs(my_rate) > target_rate) {
+                    std::cout << "Group GRAT target not met for group " << name() << std::endl;
+                    std::cout << "target = " << target_rate << '\n'
+                              << "rate   = " << my_rate << std::endl;
+                    applyProdGroupControl(ProductionSpecification::GRAT, target_rate);
+                    prodSpec().control_mode_ = ProductionSpecification::GRAT;
+                    return false;
+                }
+            }
+        }
+
+        // LRAT
+        if (prodSpec().control_mode_ != ProductionSpecification::LRAT) {
+            const double target_rate = prodSpec().liquid_max_rate_;
+            if (target_rate >= 0.0) {
+                const double my_rate = 
+                    = child_phases_summed.surf_prod_rates[phaseUsage().phase_pos[BlackoilPhases::Aqua]]
+                    + child_phases_summed.surf_prod_rates[phaseUsage().phase_pos[BlackoilPhases::Liquid]];
+                if (std::fabs(my_rate) > target_rate) {
+                    std::cout << "Group LRAT target not met for group " << name() << std::endl;
+                    std::cout << "target = " << target_rate << '\n'
+                              << "rate   = " << my_rate << std::endl;
+                    applyProdGroupControl(ProductionSpecification::LRAT, target_rate);
+                    prodSpec().control_mode_ = ProductionSpecification::LRAT;
+                    return false;
+                }
+            }
+        }
+
+        // RESV
+        if (prodSpec().control_mode_ != ProductionSpecification::RESV) {
+            const double target_rate = prodSpec().reservoir_flow_max_rate_;
+            if (target_rate >= 0.0) {
+                double my_rate = 0.0;
+                for (int phase = 0; phase < np; ++phase) {
+                    my_rate += child_phases_summed.res_prod_rates[phase];
+                }
+                if (std::fabs(my_rate) > target_rate) {
+                    std::cout << "Group RESV target not met for group " << name() << std::endl;
+                    std::cout << "target = " << target_rate << '\n'
+                              << "rate   = " << my_rate << std::endl;
+                    applyProdGroupControl(ProductionSpecification::RESV, target_rate);
+                    prodSpec().control_mode_ = ProductionSpecification::RESV;
+                    return false;
+                }
+            }
+        }
+
+
         double rate_target = std::min(std::abs(injSpec().fluid_volume_max_rate_), 
                                       prodSpec().fluid_volume_max_rate_);
         
-        double bhp_sum = child_phases_summed.bhp_sum;
         double rate_sum = child_phases_summed.rate_sum;
-        if (std::abs(bhp_sum) - std::abs(bhp_target) > epsilon) {
+        if (std::abs(rate_sum) - std::abs(rate_target) > epsilon) {
-            std::cout << "BHP not met" << std::endl;
-            std::cout << "BHP limit was " << bhp_target << std::endl;
-            std::cout << "Actual bhp was " << bhp_sum << std::endl;
-            
-            switch(prodSpec().procedure_) {
-            case ProductionSpecification::WELL:
-                getWorstOffending(well_bhp).first->shutWell();
-                return false;
-                break;
-            case ProductionSpecification::RATE:
-                applyControl(BHP);
-                return false;
-                break;
-            default:
-                // Nothing do to;
-                break;
-            }
-        }
-        if(std::abs(rate_sum) - std::abs(rate_target) > epsilon) {
             std::cout << "well_rate not met" << std::endl;
             std::cout << "target = " << rate_target 
                       << ", well_rate[index_of_well] = " 
                       << rate_sum << std::endl;
             std::cout << "Group name = " << name() << std::endl;
             
-            switch(prodSpec().procedure_) {
+            switch (prodSpec().procedure_) {
             case ProductionSpecification::WELL:
                 getWorstOffending(well_rate).first->shutWell();
                 return false;
                 break;
             case ProductionSpecification::RATE:
-                applyControl(RATE);
+                applyControl(SURFACE_RATE);
                 return false;
                 break;
             default:
@@ -231,60 +452,98 @@ namespace Opm
         }
         return max;
     }
+
+
+
+    // ==============    WellNode members   ============
+
+
     
     WellNode::WellNode(const std::string& myname,
-            ProductionSpecification prod_spec,
+                       const ProductionSpecification& prod_spec,
-            InjectionSpecification inj_spec)
+                       const InjectionSpecification& inj_spec,
-    : WellsGroupInterface(myname, prod_spec, inj_spec)
+                       const PhaseUsage& phase_usage)
+        : WellsGroupInterface(myname, prod_spec, inj_spec, phase_usage),
+          wells_(0),
+          self_index_(-1),
+          group_control_index_(-1)
     {
     }
 
-    bool WellNode::conditionsMet(const std::vector<double>& well_bhp, 
+    bool WellNode::conditionsMet(const std::vector<double>& well_bhp,
-                                 const std::vector<double>& well_rate, 
+                                 const std::vector<double>& well_reservoirrates_phase,
-                                 WellPhasesSummed& summed_phases, 
+                                 const std::vector<double>& well_surfacerates_phase,
-                                 const double epsilon)
+                                 WellPhasesSummed& summed_phases)
     {
-        
+        // Report on our rates.
-        
+        const int np = phaseUsage().num_phases;
-        // Check for self:
+        for (int phase = 0; phase < np; ++phase) {
-        if (wells_->type[self_index_] == PRODUCER) {
+            if (wells_->type[self_index_] == INJECTOR) {
-            double bhp_diff = well_bhp[self_index_] - prodSpec().BHP_limit_;
+                summed_phases.res_inj_rates[phase] = well_reservoirrates_phase[np*self_index_ + phase];
-            double rate_diff = well_rate[self_index_] - prodSpec().fluid_volume_max_rate_;
+                summed_phases.surf_inj_rates[phase] = well_surfacerates_phase[np*self_index_ + phase];
-            
+            } else {
-            if (bhp_diff > epsilon) {
+                summed_phases.res_prod_rates[phase] = well_reservoirrates_phase[np*self_index_ + phase];
-                
+                summed_phases.surf_prod_rates[phase] = well_surfacerates_phase[np*self_index_ + phase];
-                std::cout << "BHP exceeded, bhp_diff = " << bhp_diff << std::endl;
+            }
-                std::cout << "BHP_limit = " << prodSpec().BHP_limit_ << std::endl;
+        }
-                std::cout << "BHP = " << well_bhp[self_index_] << std::endl; 
+
-                shutWell();
+        // Check constraints.
-                return false;
+        bool is_producer = (wells_->type[self_index_] == PRODUCER);
+        const WellControls& ctrls = *wells_->ctrls[self_index_];
+        for (int ctrl_index = 0; ctrl_index < ctrls.num; ++ctrl_index) {
+            if (ctrl_index == ctrls.current || ctrl_index == group_control_index_) {
+                // We do not check constraints that either were used
+                // as the active control, or that come from group control.
+                continue;
+            }
+            bool ctrl_violated = false;
+            switch (ctrls.type[ctrl_index]) {
+            case BHP: {
+                const double my_well_bhp = well_bhp[self_index_];
+                const double my_target_bhp = ctrls.target[ctrl_index];
+                ctrl_violated = is_producer ? (my_target_bhp > my_well_bhp)
+                    : (my_target_bhp < my_well_bhp);
+                if (ctrl_violated) {
+                    std::cout << "BHP limit violated for well " << name() << ":\n";
+                    std::cout << "BHP limit = " << my_target_bhp << std::endl;
+                    std::cout << "BHP       = " << my_well_bhp << std::endl;
+                }
+                break;
+            }
+            case RESERVOIR_RATE: {
+                double my_rate = 0.0;
+                for (int phase = 0; phase < np; ++phase) {
+                    my_rate += ctrls.distr[np*ctrl_index + phase]*well_reservoirrates_phase[np*self_index_ + phase];
+                }
+                const double my_rate_target = ctrls.target[ctrl_index];
+                ctrl_violated = std::fabs(my_rate) > std::fabs(my_rate_target);
+                if (ctrl_violated) {
+                    std::cout << "RESERVOIR_RATE limit violated for well " << name() << ":\n";
+                    std::cout << "rate limit = " << my_rate_target;
+                    std::cout << "rate       = " << my_rate;
+                }
+                break;
+            }
+            case SURFACE_RATE: {
+                double my_rate = 0.0;
+                for (int phase = 0; phase < np; ++phase) {
+                    my_rate += ctrls.distr[np*ctrl_index + phase]*well_surfacerates_phase[np*self_index_ + phase];
+                }
+                const double my_rate_target = ctrls.target[ctrl_index];
+                ctrl_violated = std::fabs(my_rate) > std::fabs(my_rate_target);
+                if (ctrl_violated) {
+                    std::cout << "SURFACE_RATE limit violated for well " << name() << ":\n";
+                    std::cout << "rate limit = " << my_rate_target;
+                    std::cout << "rate       = " << my_rate;
+                }
+                break;
+            }
+            } // end of switch()
+            if (ctrl_violated) {
+                set_current_control(self_index_, ctrl_index, wells_);
+                return false;
             }
-            
-            if (rate_diff > epsilon) {
-                std::cout << "Rate exceeded, rate_diff = " << rate_diff << std::endl;
-                shutWell();
-                return false;
-            }
-        } else {
-            double bhp_diff = std::abs(well_bhp[self_index_]) - std::abs(injSpec().BHP_limit_);
-            double rate_diff = std::abs(well_rate[self_index_]) - std::abs(injSpec().fluid_volume_max_rate_);
-            
-            
-           if (bhp_diff > epsilon) {
-               std::cout << "BHP exceeded, bhp_diff = " << bhp_diff<<std::endl;
-               shutWell();
-               return false;
-           }
-           if (rate_diff > epsilon) {
-               std::cout << "Flow diff exceeded, flow_diff = " << rate_diff << std::endl;
-               shutWell();
-               return false;
-           }
-            
         }
-        
-        summed_phases.bhp_sum = well_bhp[self_index_];
-        summed_phases.rate_sum = well_rate[self_index_];
         return true;
     }
 
@@ -307,8 +566,14 @@ namespace Opm
     {
         wells_ = wells;
         self_index_ = self_index;
+        bool already_has_group_control = 
+            ((wells_->type[self_index_] == INJECTOR) && (injSpec().control_mode_ == InjectionSpecification::GRUP))
+            || ((wells_->type[self_index_] == PRODUCER) && (prodSpec().control_mode_ == ProductionSpecification::GRUP));
+        if (already_has_group_control) {
+            group_control_index_ = wells_->ctrls[self_index_]->num - 1;
+        }
     }
-    
+
     void WellNode::calculateGuideRates()
     {
         // Empty
@@ -328,116 +593,171 @@ namespace Opm
         return std::make_pair<WellNode*, double>(this, values[self_index_]);
     }
     
-    void WellNode::applyControl(const WellControlType type)
+    void WellNode::applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
+                                        const double target)
     {
-        wells_->ctrls[self_index_]->type[0] = type;
+        if (!wells_->type[self_index_] == INJECTOR) {
-        double target = 0.0;
+            ASSERT(target == 0.0);
-        switch(type) {
+            return;
-        case BHP:
-            if(wells_->type[self_index_] == INJECTOR) {
-                target = injSpec().BHP_limit_;
-            }
-            else {
-                target = prodSpec().BHP_limit_;
-            }
-            break;
-        case RATE:
-            if(wells_->type[self_index_] == INJECTOR) {
-                target = injSpec().fluid_volume_max_rate_;
-            }
-            else {
-                target = prodSpec().fluid_volume_max_rate_;
-            }
-            break;
         }
-        wells_->ctrls[self_index_]->target[0] = target;
+
+        const double distr[3] = { 1.0, 1.0, 1.0 };
+        WellControlType wct;
+        switch (control_mode) {
+        case InjectionSpecification::RATE:
+            wct = SURFACE_RATE;
+            break;
+        case InjectionSpecification::RESV:
+            wct = RESERVOIR_RATE;
+            break;
+        default:
+            THROW("Group injection control mode not handled: " << control_mode);
+        }
+
+        if (group_control_index_ < 0) {
+            // The well only had its own controls, no group controls.
+            append_well_controls(wct, target, distr, self_index_, wells_);
+            group_control_index_ = wells_->ctrls[self_index_]->num - 1;
+        } else {
+            // We will now modify the last control, that
+            // "belongs to" the group control.
+            const int np = wells_->number_of_phases;
+            wells_->ctrls[self_index_]->type[group_control_index_] = wct;
+            wells_->ctrls[self_index_]->target[group_control_index_] = target;
+            std::copy(distr, distr + np, wells_->ctrls[self_index_]->distr + np*group_control_index_);
+        }
+        set_current_control(self_index_, group_control_index_, wells_);
+    }
+
+
+    void WellNode::applyProdGroupControl(const ProductionSpecification::ControlMode control_mode,
+                                         const double target)
+    {
+        if (!wells_->type[self_index_] == PRODUCER) {
+            ASSERT(target == 0.0);
+            return;
+        }
+
+        double distr[3] = { 0.0, 0.0, 0.0 };
+        const int* phase_pos = phaseUsage().phase_pos;
+        const int* phase_used = phaseUsage().phase_used;
+        WellControlType wct;
+        switch (control_mode) {
+        case ProductionSpecification::ORAT:
+            wct = SURFACE_RATE;
+            if (!phase_used[BlackoilPhases::Liquid]) {
+                THROW("Oil phase not active and ORAT control specified.");
+            }
+            distr[phase_pos[BlackoilPhases::Liquid]] = 1.0;
+            break;
+        case ProductionSpecification::WRAT:
+            wct = SURFACE_RATE;
+            if (!phase_used[BlackoilPhases::Aqua]) {
+                THROW("Water phase not active and WRAT control specified.");
+            }
+            distr[phase_pos[BlackoilPhases::Aqua]] = 1.0;
+            break;
+        case ProductionSpecification::GRAT:
+            wct = SURFACE_RATE;
+            if (!phase_used[BlackoilPhases::Vapour]) {
+                THROW("Gas phase not active and GRAT control specified.");
+            }
+            distr[phase_pos[BlackoilPhases::Vapour]] = 1.0;
+            break;
+        case ProductionSpecification::LRAT:
+            wct = SURFACE_RATE;
+            if (!phase_used[BlackoilPhases::Liquid]) {
+                THROW("Oil phase not active and LRAT control specified.");
+            }
+            if (!phase_used[BlackoilPhases::Aqua]) {
+                THROW("Water phase not active and LRAT control specified.");
+            }
+            distr[phase_pos[BlackoilPhases::Liquid]] = 1.0;
+            distr[phase_pos[BlackoilPhases::Aqua]] = 1.0;
+            break;
+        case ProductionSpecification::RESV:
+            distr[0] = distr[1] = distr[2] = 1.0;
+            wct = RESERVOIR_RATE;
+            break;
+        default:
+            THROW("Group injection control mode not handled: " << control_mode);
+        }
+
+        if (group_control_index_ < 0) {
+            // The well only had its own controls, no group controls.
+            append_well_controls(wct, target, distr, self_index_, wells_);
+            group_control_index_ = wells_->ctrls[self_index_]->num - 1;
+        } else {
+            // We will now modify the last control, that
+            // "belongs to" the group control.
+            const int np = wells_->number_of_phases;
+            wells_->ctrls[self_index_]->type[group_control_index_] = wct;
+            wells_->ctrls[self_index_]->target[group_control_index_] = target;
+            std::copy(distr, distr + np, wells_->ctrls[self_index_]->distr + np*group_control_index_);
+        }
+        set_current_control(self_index_, group_control_index_, wells_);
     }
 
     namespace
     {
 
-        SurfaceComponent toSurfaceComponent(std::string type)
+        InjectionSpecification::InjectorType toInjectorType(std::string type)
         {
             if (type == "OIL") {
-                return OIL;
+                return InjectionSpecification::OIL;
             }
             if (type == "WATER") {
-                return WATER;
+                return InjectionSpecification::WATER;
             }
             if (type == "GAS") {
-                return GAS;
+                return InjectionSpecification::GAS;
             }
             THROW("Unknown type " << type << ", could not convert to SurfaceComponent");
         }
 
+
+#define HANDLE_ICM(x)                           \
+        if (type == #x) {                       \
+            return InjectionSpecification::x;   \
+        }
+
         InjectionSpecification::ControlMode toInjectionControlMode(std::string type)
         {
-            if (type == "NONE") {
+            HANDLE_ICM(NONE);
-                return InjectionSpecification::NONE;
+            HANDLE_ICM(RATE);
-            }
+            HANDLE_ICM(RESV);
+            HANDLE_ICM(BHP);
+            HANDLE_ICM(THP);
+            HANDLE_ICM(REIN);
+            HANDLE_ICM(VREP);
+            HANDLE_ICM(GRUP);
+            HANDLE_ICM(FLD);
+            THROW("Unknown type " << type << ", could not convert to InjectionSpecification::ControlMode.");
+        }
+#undef HANDLE_ICM
 
-            if (type == "ORAT") {
+#define HANDLE_PCM(x)                           \
-                return InjectionSpecification::ORAT;
+        if (type == #x) {                       \
-            }
+            return ProductionSpecification::x;  \
-            if (type == "REIN") {
-                return InjectionSpecification::REIN;
-            }
-            if (type == "RESV") {
-                return InjectionSpecification::RESV;
-            }
-            if (type == "VREP") {
-                return InjectionSpecification::VREP;
-            }
-            if (type == "WGRA") {
-                return InjectionSpecification::WGRA;
-            }
-            if (type == "FLD") {
-                return InjectionSpecification::FLD;
-            }
-            if (type == "GRUP") {
-                return InjectionSpecification::GRUP;
-            }
-
-
-            THROW("Unknown type " << type << ", could not convert to ControlMode.");
         }
 
         ProductionSpecification::ControlMode toProductionControlMode(std::string type)
         {
-            if (type == "NONE") {
+            HANDLE_PCM(NONE);
-                return ProductionSpecification::NONE_CM;
+            HANDLE_PCM(ORAT);
-            }
+            HANDLE_PCM(WRAT);
-            if (type == "ORAT") {
+            HANDLE_PCM(GRAT);
-                return ProductionSpecification::ORAT;
+            HANDLE_PCM(LRAT);
-
+            HANDLE_PCM(CRAT);
-            }
+            HANDLE_PCM(RESV);
-            if (type == "LRAT") {
+            HANDLE_PCM(PRBL);
-                return ProductionSpecification::LRAT;
+            HANDLE_PCM(BHP);
-            }
+            HANDLE_PCM(THP);
-            if (type == "REIN") {
+            HANDLE_PCM(GRUP);
-                return ProductionSpecification::REIN;
+            HANDLE_PCM(FLD);
-            }
+            THROW("Unknown type " << type << ", could not convert to ProductionSpecification::ControlMode.");
-            if (type == "RESV") {
-                return ProductionSpecification::RESV;
-            }
-            if (type == "VREP") {
-                return ProductionSpecification::VREP;
-            }
-            if (type == "WGRA") {
-                return ProductionSpecification::WGRA;
-            }
-            if (type == "FLD") {
-                return ProductionSpecification::FLD;
-            }
-            if (type == "GRUP") {
-                return ProductionSpecification::GRUP;
-            }
-            if (type == "BHP") {
-                return ProductionSpecification::BHP;
-            }
-
-            THROW("Unknown type " << type << ", could not convert to ControlMode.");
         }
+#undef HANDLE_PCM
 
         ProductionSpecification::Procedure toProductionProcedure(std::string type)
         {
@@ -456,8 +776,10 @@ namespace Opm
         }
     } // anonymous namespace
 
-    std::tr1::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name, const EclipseGridParser& deck)
+    std::tr1::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name,
+                                                               const EclipseGridParser& deck)
     {
+        PhaseUsage phase_usage = phaseUsageFromDeck(deck);
 
         std::tr1::shared_ptr<WellsGroupInterface> return_value;
         // First we need to determine whether it's a group or just a well:
@@ -482,12 +804,14 @@ namespace Opm
                     if (wconinje.wconinje[i].well_ == name) {
                         WconinjeLine line = wconinje.wconinje[i];
                         injection_specification.BHP_limit_ = line.BHP_limit_;
-                        injection_specification.injector_type_ = toSurfaceComponent(line.injector_type_);
+                        injection_specification.injector_type_ = toInjectorType(line.injector_type_);
                         injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
                         injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
-                        injection_specification.fluid_volume_max_rate_ = line.fluid_volume_max_rate_;
+                        injection_specification.reservoir_flow_max_rate_ = line.reservoir_flow_max_rate_;
                     }
                 }
+            } else {
+                injection_specification.guide_rate_ = 0.0;
             }
 
             ProductionSpecification production_specification;
@@ -498,14 +822,16 @@ namespace Opm
                     if (wconprod.wconprod[i].well_ == name) {
                         WconprodLine line = wconprod.wconprod[i];
                         production_specification.BHP_limit_ = line.BHP_limit_;
-                        production_specification.fluid_volume_max_rate_ = line.fluid_volume_max_rate_;
+                        production_specification.reservoir_flow_max_rate_ = line.reservoir_flow_max_rate_;
                         production_specification.oil_max_rate_ = line.oil_max_rate_;
                         production_specification.control_mode_ = toProductionControlMode(line.control_mode_);
-                        production_specification.water_production_target_ = line.water_max_rate_;
+                        production_specification.water_max_rate_ = line.water_max_rate_;
                     }
                 }
+            } else {
+                production_specification.guide_rate_ = 0.0;
             }
-            return_value.reset(new WellNode(name, production_specification, injection_specification));
+            return_value.reset(new WellNode(name, production_specification, injection_specification, phase_usage));
         } else {
             InjectionSpecification injection_specification;
             if (deck.hasField("GCONINJE")) {
@@ -513,10 +839,10 @@ namespace Opm
                 for (size_t i = 0; i < gconinje.gconinje.size(); i++) {
                     if (gconinje.gconinje[i].group_ == name) {
                         GconinjeLine line = gconinje.gconinje[i];
-                        injection_specification.injector_type_ = toSurfaceComponent(line.injector_type_);
+                        injection_specification.injector_type_ = toInjectorType(line.injector_type_);
                         injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_);
                         injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_;
-                        injection_specification.fluid_volume_max_rate_ = line.resv_flow_max_rate_;
+                        injection_specification.reservoir_flow_max_rate_ = line.resv_flow_max_rate_;
                     }
                 }
             }
@@ -532,7 +858,7 @@ namespace Opm
                         production_specification.oil_max_rate_ = line.oil_max_rate_;
                         std::cout << "control_mode = " << line.control_mode_ << std::endl;
                         production_specification.control_mode_ = toProductionControlMode(line.control_mode_);
-                        production_specification.water_production_target_ = line.water_max_rate_;
+                        production_specification.water_max_rate_ = line.water_max_rate_;
                         production_specification.gas_max_rate_ = line.gas_max_rate_;
                         production_specification.liquid_max_rate_ = line.liquid_max_rate_;
                         production_specification.procedure_ = toProductionProcedure(line.procedure_);
@@ -540,7 +866,7 @@ namespace Opm
                 }
             }
 
-            return_value.reset(new WellsGroup(name, production_specification, injection_specification));
+            return_value.reset(new WellsGroup(name, production_specification, injection_specification, phase_usage));
         }
 
         return return_value;

opm/core/WellsGroup.hpp

@@ -1,3 +1,22 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  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_WELLSGROUP_HPP
 #define	OPM_WELLSGROUP_HPP
 
@@ -5,6 +24,7 @@
 #include <opm/core/ProductionSpecification.hpp>
 #include <opm/core/eclipse/EclipseGridParser.hpp>
 #include <opm/core/grid.h>
+#include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
 #include <string>
 
 
@@ -16,11 +36,14 @@ namespace Opm
     
     /// Basic information needed for group control (each group should typically
     /// not exceed the sum of its leaf nodes)
-    struct WellPhasesSummed {
+    struct WellPhasesSummed
+    {
         WellPhasesSummed();
-        double bhp_sum;
+        double res_inj_rates[3];
-        double rate_sum;
+        double res_prod_rates[3];
-        
+        double surf_inj_rates[3];
+        double surf_prod_rates[3];
+
         /// Sums each component
         void operator+=(const WellPhasesSummed& other);
     };
@@ -29,8 +52,9 @@ namespace Opm
     {
     public:
         WellsGroupInterface(const std::string& name,
-                            ProductionSpecification prod_spec,
+                            const ProductionSpecification& prod_spec,
-                            InjectionSpecification inj_spec);
+                            const InjectionSpecification& inj_spec,
+                            const PhaseUsage& phase_usage);
         virtual ~WellsGroupInterface();
 
         /// The unique identifier for the well or well group.
@@ -47,7 +71,9 @@ namespace Opm
         
         /// Injection specifications for the well or well group.
         InjectionSpecification& injSpec();
-        
+
+        /// Phase usage information.
+        const PhaseUsage& phaseUsage() const;
         
         /// \returns true if the object is a leaf node (WellNode), false otherwise.
         virtual bool isLeafNode() const;
@@ -76,7 +102,7 @@ namespace Opm
         /// down wells). Only one change is applied per invocation. Typical use will be
         /// \code
         /// solve_pressure();
-        /// while(!group.conditionsMet(well_bhp, well_rate, summed_phases)) {
+        /// while(!group.conditionsMet(...)) {
         ///     solve_pressure();
         /// }
         /// \endcode
@@ -84,23 +110,33 @@ namespace Opm
         /// \note It's highly recommended to use the conditionsMet found in WellsManager.
         /// \param[in]    well_bhp  A vector containing the bhp for each well. Is assumed 
         ///                         to be ordered the same way as the related Wells-struct.
-        /// \param[in]    well_rate A vector containing the rate for each well. Is assumed 
+        /// \param[in]    well_reservoirrates_phase
-        ///                         to be ordered the same way as the related Wells-struct.
+        ///                         A vector containing reservoir rates by phase for each well.
-        /// \param[out]   summed_phases Will at end of invocation contain the summed phases
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
-        ///                             (bhp, rate ,etc.) for the group.
+        ///                         with all phase rates of a single well adjacent in the array.
-        /// \param[in]    epsilon   The error tolerated for each inequality. Formally, it will accept
+        /// \param[in]    well_surfacerates_phase
-        ///                         (a - b <= epsilon) as (a <= b).
+        ///                         A vector containing surface rates by phase for each well.
+        ///                         Is assumed to be ordered the same way as the related Wells-struct,
+        ///                         with all phase rates of a single well adjacent in the array.
+        /// \param[out]   summed_phases Will at end of invocation contain the summed phase rates
+        ///                             (rate ,etc.) for the group.
         /// \return true if no violations were found, false otherwise (false also implies a change).
         virtual bool conditionsMet(const std::vector<double>& well_bhp,
-                                   const std::vector<double>& well_rate,
+                                   const std::vector<double>& well_reservoirrates_phase,
-                                   WellPhasesSummed& summed_phases,
+                                   const std::vector<double>& well_surfacerates_phase,
-                                   const double epsilon = 1e-8) = 0;
+                                   WellPhasesSummed& summed_phases);
-        
-        /// Sets the current active control to the provided one for all wells within the group
-        /// \note Also changes the target based on type.
-        /// \param[in] type the type to change to which the control is changed.
-        virtual void applyControl(const WellControlType type) = 0;
         
+        /// Sets the current active control to the provided one for all injectors within the group.
+        /// After this call, the combined rate (which rate depending on control_mode) of the group
+        /// shall be equal to target.
+        virtual void applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
+                                          const double target) = 0;
+        /// Sets the current active control to the provided one for all producers within the group.
+        /// After this call, the combined rate (which rate depending on control_mode) of the group
+        /// shall be equal to target.
+        virtual void applyProdGroupControl(const ProductionSpecification::ControlMode control_mode,
+                                           const double target) = 0;
+
         /// Gets the worst offending well based on the input
         /// \param values A vector of a values for each well. This is assumed to be ordered the same way as the 
         ///               relevant Wells struct.
@@ -108,12 +144,23 @@ namespace Opm
         virtual std::pair<WellNode*, double> getWorstOffending(const std::vector<double>& values) = 0;
         
     protected:
-           WellsGroupInterface* parent_;
+        /// Calculates the correct rate for the given ProductionSpecification::ControlMode
+        double rateByMode(const double* res_rates, 
+                          const double* surf_rates,
+                          const ProductionSpecification::ControlMode mode);
+
+        /// Calculates the correct rate for the given InjectionSpecification::ControlMode
+        double rateByMode(const double* res_rates, 
+                          const double* surf_rates,
+                          const InjectionSpecification::ControlMode mode);
+
+        WellsGroupInterface* parent_;
 
     private:
         std::string name_;
         ProductionSpecification production_specification_;
         InjectionSpecification injection_specification_;
+        PhaseUsage phase_usage_;
     };
 
 
@@ -122,24 +169,36 @@ namespace Opm
     {
     public:
         WellsGroup(const std::string& name,
-                   ProductionSpecification prod_spec,
+                   const ProductionSpecification& prod_spec,
-                   InjectionSpecification inj_spec);
+                   const InjectionSpecification& inj_spec,
+                   const PhaseUsage& phase_usage);
 
         virtual WellsGroupInterface* findGroup(const std::string& name_of_node);
 
         void addChild(std::tr1::shared_ptr<WellsGroupInterface> child);
         
         virtual bool conditionsMet(const std::vector<double>& well_bhp,
-                                   const std::vector<double>& well_rate,
+                                   const std::vector<double>& well_reservoirrates_phase,
-                                   WellPhasesSummed& summed_phases,
+                                   const std::vector<double>& well_surfacerates_phase,
-                                   const double epsilon = 1e-8);
+                                   WellPhasesSummed& summed_phases);
         
         
         virtual void calculateGuideRates();
 
         virtual int numberOfLeafNodes();
         virtual std::pair<WellNode*, double> getWorstOffending(const std::vector<double>& values);
-        virtual void applyControl(const WellControlType type);
+
+        /// Sets the current active control to the provided one for all injectors within the group.
+        /// After this call, the combined rate (which rate depending on control_mode) of the group
+        /// shall be equal to target.
+        virtual void applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
+                                          const double target);
+
+        /// Sets the current active control to the provided one for all producers within the group.
+        /// After this call, the combined rate (which rate depending on control_mode) of the group
+        /// shall be equal to target.
+        virtual void applyProdGroupControl(const ProductionSpecification::ControlMode control_mode,
+                                           const double target);
 
     private:
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > children_;
@@ -151,14 +210,15 @@ namespace Opm
     {
     public:
         WellNode(const std::string& name,
-                 ProductionSpecification prod_spec,
+                 const ProductionSpecification& prod_spec,
-                InjectionSpecification inj_spec);
+                 const InjectionSpecification& inj_spec,
+                 const PhaseUsage& phase_usage);
 
         virtual WellsGroupInterface* findGroup(const std::string& name_of_node);
         virtual bool conditionsMet(const std::vector<double>& well_bhp,
-                                   const std::vector<double>& well_rate,
+                                   const std::vector<double>& well_reservoirrates_phase,
-                                   WellPhasesSummed& summed_phases,
+                                   const std::vector<double>& well_surfacerates_phase,
-                                   const double epsilon = 1e-8);
+                                   WellPhasesSummed& summed_phases);
         
         virtual bool isLeafNode() const;
         
@@ -171,11 +231,23 @@ namespace Opm
         void shutWell();
         
         virtual std::pair<WellNode*, double> getWorstOffending(const std::vector<double>& values);
-        virtual void applyControl(const WellControlType type);
+
-        
+        /// Sets the current active control to the provided one for all injectors within the group.
+        /// After this call, the combined rate (which rate depending on control_mode) of the group
+        /// shall be equal to target.
+        virtual void applyInjGroupControl(const InjectionSpecification::ControlMode control_mode,
+                                          const double target);
+
+        /// Sets the current active control to the provided one for all producers within the group.
+        /// After this call, the combined rate (which rate depending on control_mode) of the group
+        /// shall be equal to target.
+        virtual void applyProdGroupControl(const ProductionSpecification::ControlMode control_mode,
+                                           const double target);
+
     private:
         Wells* wells_;
         int self_index_;
+        int group_control_index_;
     };
 
     /// Creates the WellsGroupInterface for the given name