diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index d155afb87..471bab68b 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -1,36 +1,64 @@ -/* - * File: WellsGroup.cpp - * Author: kjetilo - * - * Created on March 27, 2012, 9:27 AM - */ +/* + 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 . +*/ #include #include #include +#include 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; - bhp_sum += other.bhp_sum; + for (int i = 0; i < 3; ++i) { + 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, - InjectionSpecification inje_spec) - : parent_(NULL), - name_(myname), - production_specification_(prod_spec), - injection_specification_(inje_spec) + const ProductionSpecification& prod_spec, + const InjectionSpecification& inje_spec, + const PhaseUsage& phase_usage) + : parent_(NULL), + name_(myname), + 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& well_bhp, - const std::vector& well_rate, - WellPhasesSummed& summed_phases, - const double epsilon) + + void WellsGroup::calculateGuideRates() { + 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& well_bhp, + const std::vector& well_reservoirrates_phase, + const std::vector& 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, - current_child_phases_summed, epsilon)) { + if (!children_[i]->conditionsMet(well_bhp, + 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) { - 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) { + 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, - InjectionSpecification inj_spec) - : WellsGroupInterface(myname, prod_spec, inj_spec) + const ProductionSpecification& prod_spec, + const InjectionSpecification& 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& well_bhp, - const std::vector& well_rate, - WellPhasesSummed& summed_phases, - const double epsilon) + bool WellNode::conditionsMet(const std::vector& well_bhp, + const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase, + WellPhasesSummed& summed_phases) { - - - // Check for self: - if (wells_->type[self_index_] == PRODUCER) { - double bhp_diff = well_bhp[self_index_] - prodSpec().BHP_limit_; - double rate_diff = well_rate[self_index_] - prodSpec().fluid_volume_max_rate_; - - if (bhp_diff > epsilon) { - - 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(); - return false; + // Report on our rates. + const int np = phaseUsage().num_phases; + for (int phase = 0; phase < np; ++phase) { + if (wells_->type[self_index_] == INJECTOR) { + summed_phases.res_inj_rates[phase] = well_reservoirrates_phase[np*self_index_ + phase]; + summed_phases.surf_inj_rates[phase] = well_surfacerates_phase[np*self_index_ + phase]; + } else { + 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]; + } + } + + // Check constraints. + 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< 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(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; - double target = 0.0; - switch(type) { - 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; + if (!wells_->type[self_index_] == INJECTOR) { + ASSERT(target == 0.0); + return; } - 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") { - return InjectionSpecification::NONE; - } + HANDLE_ICM(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") { - return InjectionSpecification::ORAT; - } - 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."); +#define HANDLE_PCM(x) \ + if (type == #x) { \ + return ProductionSpecification::x; \ } ProductionSpecification::ControlMode toProductionControlMode(std::string type) { - if (type == "NONE") { - return ProductionSpecification::NONE_CM; - } - if (type == "ORAT") { - return ProductionSpecification::ORAT; - - } - if (type == "LRAT") { - return ProductionSpecification::LRAT; - } - if (type == "REIN") { - return ProductionSpecification::REIN; - } - 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."); + HANDLE_PCM(NONE); + HANDLE_PCM(ORAT); + HANDLE_PCM(WRAT); + HANDLE_PCM(GRAT); + HANDLE_PCM(LRAT); + HANDLE_PCM(CRAT); + HANDLE_PCM(RESV); + HANDLE_PCM(PRBL); + HANDLE_PCM(BHP); + HANDLE_PCM(THP); + HANDLE_PCM(GRUP); + HANDLE_PCM(FLD); + THROW("Unknown type " << type << ", could not convert to ProductionSpecification::ControlMode."); } +#undef HANDLE_PCM ProductionSpecification::Procedure toProductionProcedure(std::string type) { @@ -456,8 +776,10 @@ namespace Opm } } // anonymous namespace - std::tr1::shared_ptr createWellsGroup(const std::string& name, const EclipseGridParser& deck) + std::tr1::shared_ptr createWellsGroup(const std::string& name, + const EclipseGridParser& deck) { + PhaseUsage phase_usage = phaseUsageFromDeck(deck); std::tr1::shared_ptr 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; diff --git a/opm/core/WellsGroup.hpp b/opm/core/WellsGroup.hpp index f31f3aff9..f0f7c7901 100644 --- a/opm/core/WellsGroup.hpp +++ b/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 . +*/ + #ifndef OPM_WELLSGROUP_HPP #define OPM_WELLSGROUP_HPP @@ -5,6 +24,7 @@ #include #include #include +#include #include @@ -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 rate_sum; - + double res_inj_rates[3]; + 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, - InjectionSpecification inj_spec); + const ProductionSpecification& prod_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 - /// to be ordered the same way as the related Wells-struct. - /// \param[out] summed_phases Will at end of invocation contain the summed phases - /// (bhp, rate ,etc.) for the group. - /// \param[in] epsilon The error tolerated for each inequality. Formally, it will accept - /// (a - b <= epsilon) as (a <= b). + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir 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[in] well_surfacerates_phase + /// 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& well_bhp, - const std::vector& well_rate, - WellPhasesSummed& summed_phases, - const double epsilon = 1e-8) = 0; - - /// 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; + const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase, + WellPhasesSummed& summed_phases); + /// 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 getWorstOffending(const std::vector& 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, - InjectionSpecification inj_spec); + const ProductionSpecification& prod_spec, + const InjectionSpecification& inj_spec, + const PhaseUsage& phase_usage); virtual WellsGroupInterface* findGroup(const std::string& name_of_node); void addChild(std::tr1::shared_ptr child); virtual bool conditionsMet(const std::vector& well_bhp, - const std::vector& well_rate, - WellPhasesSummed& summed_phases, - const double epsilon = 1e-8); + const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase, + WellPhasesSummed& summed_phases); virtual void calculateGuideRates(); virtual int numberOfLeafNodes(); virtual std::pair getWorstOffending(const std::vector& 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 > children_; @@ -151,14 +210,15 @@ namespace Opm { public: WellNode(const std::string& name, - ProductionSpecification prod_spec, - InjectionSpecification inj_spec); + const ProductionSpecification& prod_spec, + const InjectionSpecification& inj_spec, + const PhaseUsage& phase_usage); virtual WellsGroupInterface* findGroup(const std::string& name_of_node); virtual bool conditionsMet(const std::vector& well_bhp, - const std::vector& well_rate, - WellPhasesSummed& summed_phases, - const double epsilon = 1e-8); + const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase, + WellPhasesSummed& summed_phases); virtual bool isLeafNode() const; @@ -171,11 +231,23 @@ namespace Opm void shutWell(); virtual std::pair getWorstOffending(const std::vector& 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