OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3334 from akva2/wellif_further_split

Browse Source 
Further splitting in WellInterface
This commit is contained in:
Bård Skaflestad 2021-06-07 21:41:12 +02:00 committed by GitHub
commit 4767311f76
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
14 changed files with 1035 additions and 626 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -66,8 +66,10 @@ list (APPEND MAIN_SOURCE_FILES
opm/simulators/wells/VFPProdProperties.cpp
opm/simulators/wells/VFPInjProperties.cpp
opm/simulators/wells/WellGroupHelpers.cpp
opm/simulators/wells/WellInterfaceEval.cpp
opm/simulators/wells/WellInterfaceFluidSystem.cpp
opm/simulators/wells/WellInterfaceGeneric.cpp
opm/simulators/wells/WellInterfaceIndices.cpp
opm/simulators/wells/WellProdIndexCalculator.cpp
opm/simulators/wells/WellState.cpp
opm/simulators/wells/WGState.cpp

2
opm/simulators/wells/BlackoilWellModel_impl.hpp
View File

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

4
opm/simulators/wells/MultisegmentWell_impl.hpp
View File

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

1
opm/simulators/wells/StandardWell.hpp
View File

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

4
opm/simulators/wells/StandardWell_impl.hpp
View File

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

74
opm/simulators/wells/WellInterface.hpp
View File

@ -30,13 +30,9 @@
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTestState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellGroupHelpers.hpp>
#include <opm/simulators/wells/WellProdIndexCalculator.hpp>
#include <opm/simulators/wells/WellState.hpp>
// NOTE: GasLiftSingleWell.hpp includes StandardWell.hpp which includes ourself
@ -56,22 +52,20 @@ namespace Opm {
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/matrixmatrix.hh>
#include <opm/material/densead/Math.hpp>
#include <opm/material/densead/Evaluation.hpp>
#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <array>
#include <cassert>
#include <memory>
#include <string>
#include <vector>
namespace Opm
{
template<typename TypeTag>
class WellInterface : public WellInterfaceFluidSystem<GetPropType<TypeTag, Properties::FluidSystem>>
class WellInterface : public WellInterfaceIndices<GetPropType<TypeTag, Properties::FluidSystem>,
GetPropType<TypeTag, Properties::Indices>,
GetPropType<TypeTag, Properties::Scalar>>
{
public:
@ -230,7 +224,6 @@ public:
Scalar volumetricSurfaceRateForConnection(int cellIdx, int phaseIdx) const;
template <class EvalWell>
Eval restrictEval(const EvalWell& in) const
{
@ -292,28 +285,17 @@ protected:
bool changed_to_stopped_this_step_ = false;
int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const;
int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
double wpolymer() const;
double wfoam() const;
double wsalt() const;
template <class ValueType>
ValueType calculateBhpFromThp(const WellState& well_state, const std::vector<ValueType>& rates, const Well& well, const SummaryState& summaryState, DeferredLogger& deferred_logger) const;
virtual double getRefDensity() const = 0;
// Component fractions for each phase for the well
const std::vector<double>& compFrac() const;
double scalingFactor(const int comp_idx) const;
std::vector<double> initialWellRateFractions(const Simulator& ebosSimulator, const WellState& well_state) const;
// check whether the well is operable under BHP limit with current reservoir condition
@ -361,54 +343,6 @@ protected:
void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
DeferredLogger& deferred_logger);
template <class EvalWell>
void getGroupInjectionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
const EvalWell& bhp,
const EvalWell& injection_rate,
EvalWell& control_eq,
double efficiencyFactor,
DeferredLogger& deferred_logger);
template <class EvalWell>
void getGroupProductionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const EvalWell& bhp,
const std::vector<EvalWell>& rates,
EvalWell& control_eq,
double efficiencyFactor);
template <class EvalWell, class BhpFromThpFunc>
void assembleControlEqInj(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::InjectionControls& controls,
const EvalWell& bhp,
const EvalWell& injection_rate,
BhpFromThpFunc bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger);
template <class EvalWell, class BhpFromThpFunc>
void assembleControlEqProd(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::ProductionControls& controls,
const EvalWell& bhp,
const std::vector<EvalWell>& rates, // Always 3 canonical rates.
BhpFromThpFunc bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger);
};
}

572
opm/simulators/wells/WellInterfaceEval.cpp Normal file
View File

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

166
opm/simulators/wells/WellInterfaceEval.hpp Normal file
View File

@ -0,0 +1,166 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2017 IRIS
Copyright 2019 Norce
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLINTERFACE_EVAL_HEADER_INCLUDED
#define OPM_WELLINTERFACE_EVAL_HEADER_INCLUDED
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/ScheduleTypes.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <functional>
namespace Opm
{
class DeferredLogger;
class Group;
class GroupState;
class Schedule;
class SummaryState;
template<class FluidSystem> class WellInterfaceFluidSystem;
class WellState;
template<class FluidSystem>
class WellInterfaceEval {
static constexpr int Water = BlackoilPhases::Aqua;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Gas = BlackoilPhases::Vapour;
protected:
WellInterfaceEval(const WellInterfaceFluidSystem<FluidSystem>& baseif);
template<class EvalWell>
void getGroupInjectionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const InjectorType& injectorType,
const EvalWell& bhp,
const EvalWell& injection_rate,
EvalWell& control_eq,
double efficiencyFactor,
DeferredLogger& deferred_logger);
template<class EvalWell>
void getGroupProductionControl(const Group& group,
const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const EvalWell& bhp,
const std::vector<EvalWell>& rates,
EvalWell& control_eq,
double efficiencyFactor);
template<class EvalWell, class BhpFromThpFunc>
void assembleControlEqProd(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::ProductionControls& controls,
const EvalWell& bhp,
const std::vector<EvalWell>& rates, // Always 3 canonical rates.
BhpFromThpFunc bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger)
{
std::function<EvalWell()> eval = [&bhp_from_thp]() { return bhp_from_thp(); };
assembleControlEqProd_(well_state,
group_state,
schedule,
summaryState,
controls,
bhp,
rates,
eval,
control_eq,
deferred_logger);
}
template<class EvalWell>
void assembleControlEqProd_(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::ProductionControls& controls,
const EvalWell& bhp,
const std::vector<EvalWell>& rates, // Always 3 canonical rates.
const std::function<EvalWell()>& bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger);
template<class EvalWell, class BhpFromThpFunc>
void assembleControlEqInj(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::InjectionControls& controls,
const EvalWell& bhp,
const EvalWell& injection_rate,
BhpFromThpFunc bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger)
{
std::function<EvalWell()> eval = [&bhp_from_thp]() { return bhp_from_thp(); };
assembleControlEqInj_(well_state,
group_state,
schedule,
summaryState,
controls,
bhp,
injection_rate,
eval,
control_eq,
deferred_logger);
}
template<class EvalWell>
void assembleControlEqInj_(const WellState& well_state,
const GroupState& group_state,
const Schedule& schedule,
const SummaryState& summaryState,
const Well::InjectionControls& controls,
const EvalWell& bhp,
const EvalWell& injection_rate,
const std::function<EvalWell()>& bhp_from_thp,
EvalWell& control_eq,
DeferredLogger& deferred_logger);
template <class EvalWell>
EvalWell calculateBhpFromThp(const WellState& well_state,
const std::vector<EvalWell>& rates,
const Well& well,
const SummaryState& summaryState,
const double rho,
DeferredLogger& deferred_logger) const;
const WellInterfaceFluidSystem<FluidSystem>& baseif_;
};
}
#endif // OPM_WELLINTERFACE_EVAL_HEADER_INCLUDED

17
opm/simulators/wells/WellInterfaceFluidSystem.cpp
View File

@ -899,6 +899,23 @@ checkMaxRatioLimitWell(const WellState& well_state,
return (well_ratio > max_ratio_limit);
}
template<typename FluidSystem>
int
WellInterfaceFluidSystem<FluidSystem>::
flowPhaseToEbosPhaseIdx(const int phaseIdx) const
{
const auto& pu = this->phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx)
return FluidSystem::waterPhaseIdx;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx)
return FluidSystem::oilPhaseIdx;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx)
return FluidSystem::gasPhaseIdx;
// for other phases return the index
return phaseIdx;
}
template class WellInterfaceFluidSystem<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>>;
template class WellInterfaceFluidSystem<BlackOilFluidSystem<double,EclAlternativeBlackOilIndexTraits>>;

17
opm/simulators/wells/WellInterfaceFluidSystem.hpp
View File

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

33
opm/simulators/wells/WellInterfaceGeneric.hpp
View File

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

147
opm/simulators/wells/WellInterfaceIndices.cpp Normal file
View File

@ -0,0 +1,147 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2018 IRIS
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#include <opm/simulators/wells/WellInterfaceIndices.hpp>
#include <ebos/eclalternativeblackoilindices.hh>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/models/blackoil/blackoilindices.hh>
#include <opm/models/blackoil/blackoilonephaseindices.hh>
#include <opm/models/blackoil/blackoiltwophaseindices.hh>
#include <cassert>
namespace Opm
{
template<class FluidSystem, class Indices, class Scalar>
WellInterfaceIndices<FluidSystem,Indices,Scalar>::
WellInterfaceIndices(const Well& well,
const ParallelWellInfo& parallel_well_info,
const int time_step,
const typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data)
: WellInterfaceFluidSystem<FluidSystem>(well,
parallel_well_info,
time_step,
rate_converter,
pvtRegionIdx,
num_components,
num_phases,
index_of_well,
first_perf_index,
perf_data)
, WellInterfaceEval<FluidSystem>(static_cast<const WellInterfaceFluidSystem<FluidSystem>&>(*this))
{
}
template<class FluidSystem, class Indices, class Scalar>
int
WellInterfaceIndices<FluidSystem,Indices,Scalar>::
flowPhaseToEbosCompIdx(const int phaseIdx) const
{
const auto& pu = this->phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx)
return Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx)
return Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx)
return Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
// for other phases return the index
return phaseIdx;
}
template<class FluidSystem, class Indices, class Scalar>
int
WellInterfaceIndices<FluidSystem,Indices,Scalar>::
ebosCompIdxToFlowCompIdx(const unsigned compIdx) const
{
const auto& pu = this->phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx) == compIdx)
return pu.phase_pos[Water];
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx) == compIdx)
return pu.phase_pos[Oil];
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx) == compIdx)
return pu.phase_pos[Gas];
// for other phases return the index
return compIdx;
}
template<class FluidSystem, class Indices, class Scalar>
double
WellInterfaceIndices<FluidSystem,Indices,Scalar>::
scalingFactor(const int phaseIdx) const
{
const auto& pu = this->phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx)
return 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && pu.phase_pos[Oil] == phaseIdx)
return 1.0;
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && pu.phase_pos[Gas] == phaseIdx)
return 0.01;
if (Indices::enableSolvent && phaseIdx == Indices::contiSolventEqIdx )
return 0.01;
// we should not come this far
assert(false);
return 1.0;
}
#define INSTANCE(A, ...) \
template class WellInterfaceIndices<BlackOilFluidSystem<double, A>, \
__VA_ARGS__, \
double>;
// One phase
INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilOnePhaseIndices<0u,0u,0u,1u,false,false,0u,1u>)
// Two phase
INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,0u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,1u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,false,0u,2u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,1u,0u,false,false,0u,2u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,2u,0u,false,false,0u,2u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilTwoPhaseIndices<0u,0u,0u,0u,false,true,0u,2u>)
// Blackoil
INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,true,false,0u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,true,0u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<1u,0u,0u,0u,false,false,0u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,1u,0u,0u,false,false,0u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,1u,0u,false,false,0u>)
INSTANCE(BlackOilDefaultIndexTraits,BlackOilIndices<0u,0u,0u,1u,false,false,0u>)
// Alternative indices
INSTANCE(EclAlternativeBlackOilIndexTraits,BlackOilIndices<0u,0u,0u,0u,false,false,0u>)
} // namespace Opm

61
opm/simulators/wells/WellInterfaceIndices.hpp Normal file
View File

@ -0,0 +1,61 @@
/*
Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
Copyright 2017 Statoil ASA.
Copyright 2017 IRIS
Copyright 2019 Norce
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WELLINTERFACE_INDICES_HEADER_INCLUDED
#define OPM_WELLINTERFACE_INDICES_HEADER_INCLUDED
#include <opm/simulators/wells/WellInterfaceFluidSystem.hpp>
#include <opm/simulators/wells/WellInterfaceEval.hpp>
namespace Opm
{
template<class FluidSystem, class Indices, class Scalar>
class WellInterfaceIndices : public WellInterfaceFluidSystem<FluidSystem>
, public WellInterfaceEval<FluidSystem>
{
public:
using WellInterfaceFluidSystem<FluidSystem>::Gas;
using WellInterfaceFluidSystem<FluidSystem>::Oil;
using WellInterfaceFluidSystem<FluidSystem>::Water;
protected:
WellInterfaceIndices(const Well& well,
const ParallelWellInfo& parallel_well_info,
const int time_step,
const typename WellInterfaceFluidSystem<FluidSystem>::RateConverterType& rate_converter,
const int pvtRegionIdx,
const int num_components,
const int num_phases,
const int index_of_well,
const int first_perf_index,
const std::vector<PerforationData>& perf_data);
int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
double scalingFactor(const int phaseIdx) const;
};
}
#endif // OPM_WELLINTERFACE_INDICES_HEADER_INCLUDED

561
opm/simulators/wells/WellInterface_impl.hpp
View File

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