OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2026-09-05 04:40:19 -05:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #3487 from joakim-hove/single-well-state

Browse Source 
Single well state
This commit is contained in:
Joakim Hove
2021-08-25 14:27:41 +02:00
committed by GitHub
parent 36201c8871 94d12e1230
commit 1a02ac7bd5
22 changed files with 363 additions and 485 deletions
Download Patch File Download Diff File Expand all files Collapse all files
ebos/ecltracermodel.hh
+2 -2
View File
@@ -433,8 +433,8 @@ protected:
// TODO: Some inconsistencies here that perhaps should be clarified. The "offical" rate as reported below is
// occasionally significant different from the sum over connections (as calculated above). Only observed
// for small values, neglible for the rate itself, but matters when used to calculate tracer concentrations.
std::size_t well_index = simulator_.problem().wellModel().wellState().wellIndex(well.name());
Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().wellRates(well_index)[tr.phaseIdx_];
std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value();
Scalar official_well_rate_total = simulator_.problem().wellModel().wellState().well(well_index).surface_rates[tr.phaseIdx_];
rateWellTotal = official_well_rate_total;
opm/simulators/wells/BlackoilWellModelGeneric.cpp
+21 -22
View File
@@ -169,9 +169,9 @@ loadRestartData(const data::Wells& rst_wells,
phs.at( phases.phase_pos[BlackoilPhases::Vapour] ) = rt::gas;
}
for( const auto& wm : well_state.wellMap() ) {
const auto well_index = wm.second[ 0 ];
const auto& rst_well = rst_wells.at( wm.first );
for( std::size_t well_index = 0; well_index < well_state.size(); well_index++) {
const auto& well_name = well_state.name(well_index);
const auto& rst_well = rst_wells.at(well_name);
auto& ws = well_state.well(well_index);
ws.bhp = rst_well.bhp;
ws.thp = rst_well.thp;
@@ -186,10 +186,10 @@ loadRestartData(const data::Wells& rst_wells,
for( size_t i = 0; i < phs.size(); ++i ) {
assert( rst_well.rates.has( phs[ i ] ) );
well_state.wellRates(well_index)[i] = rst_well.rates.get(phs[i]);
ws.surface_rates[i] = rst_well.rates.get(phs[i]);
}
auto& perf_data = well_state.perfData(well_index);
auto& perf_data = ws.perf_data;
auto& perf_pressure = perf_data.pressure;
auto& perf_rates = perf_data.rates;
auto& perf_phase_rates = perf_data.phase_rates;
@@ -206,7 +206,6 @@ loadRestartData(const data::Wells& rst_wells,
if (handle_ms_well && !rst_well.segments.empty()) {
// we need the well_ecl_ information
const std::string& well_name = wm.first;
const Well& well_ecl = getWellEcl(well_name);
const WellSegments& segment_set = well_ecl.getSegments();
@@ -447,7 +446,7 @@ checkGroupInjectionConstraints(const Group& group,
if (currentControl != Group::InjectionCMode::RATE)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
current_rate += WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
// sum over all nodes
current_rate = comm_.sum(current_rate);
@@ -483,13 +482,13 @@ checkGroupInjectionConstraints(const Group& group,
{
double production_Rate = 0.0;
const Group& groupRein = schedule().getGroup(controls.reinj_group, reportStepIdx);
production_Rate += WellGroupHelpers::sumWellRates(groupRein, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false);
production_Rate += WellGroupHelpers::sumWellSurfaceRates(groupRein, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/false);
// sum over all nodes
production_Rate = comm_.sum(production_Rate);
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
current_rate += WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, phasePos, /*isInjector*/true);
// sum over all nodes
current_rate = comm_.sum(current_rate);
@@ -550,7 +549,7 @@ checkGroupProductionConstraints(const Group& group,
if (currentControl != Group::ProductionCMode::ORAT)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
current_rate += WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
// sum over all nodes
current_rate = comm_.sum(current_rate);
@@ -570,7 +569,7 @@ checkGroupProductionConstraints(const Group& group,
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
current_rate += WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
// sum over all nodes
current_rate = comm_.sum(current_rate);
@@ -588,7 +587,7 @@ checkGroupProductionConstraints(const Group& group,
if (currentControl != Group::ProductionCMode::GRAT)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
current_rate += WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Vapour], false);
// sum over all nodes
current_rate = comm_.sum(current_rate);
@@ -605,8 +604,8 @@ checkGroupProductionConstraints(const Group& group,
if (currentControl != Group::ProductionCMode::LRAT)
{
double current_rate = 0.0;
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
current_rate += WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
current_rate += WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Liquid], false);
current_rate += WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, phase_usage_.phase_pos[BlackoilPhases::Aqua], false);
// sum over all nodes
current_rate = comm_.sum(current_rate);
@@ -678,8 +677,8 @@ checkGconsaleLimits(const Group& group,
const Group::ProductionCMode& oldProductionControl = this->groupState().production_control(group.name());
int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
double production_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/false);
double injection_rate = WellGroupHelpers::sumWellRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/true);
double production_rate = WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/false);
double injection_rate = WellGroupHelpers::sumWellSurfaceRates(group, schedule(), well_state, reportStepIdx, gasPos, /*isInjector*/true);
// sum over all nodes
injection_rate = comm_.sum(injection_rate);
@@ -807,7 +806,7 @@ checkGroupHigherConstraints(const Group& group,
calcInjRates(fipnum, pvtreg, resv_coeff_inj);
for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) {
const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ true);
const double local_current_rate = WellGroupHelpers::sumWellSurfaceRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ true);
// Sum over all processes
rates[phasePos] = comm_.sum(local_current_rate);
}
@@ -845,7 +844,7 @@ checkGroupHigherConstraints(const Group& group,
if (!skip && group.isProductionGroup()) {
// Obtain rates for group.
for (int phasePos = 0; phasePos < phase_usage_.num_phases; ++phasePos) {
const double local_current_rate = WellGroupHelpers::sumWellRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ false);
const double local_current_rate = WellGroupHelpers::sumWellSurfaceRates(group, schedule(), this->wellState(), reportStepIdx, phasePos, /* isInjector */ false);
// Sum over all processes
rates[phasePos] = -comm_.sum(local_current_rate);
}
@@ -1132,7 +1131,7 @@ updateWsolvent(const Group& group,
int gasPos = phase_usage_.phase_pos[BlackoilPhases::Vapour];
const auto& controls = group.injectionControls(Phase::GAS, summaryState_);
const Group& groupRein = schedule_.getGroup(controls.reinj_group, reportStepIdx);
double gasProductionRate = WellGroupHelpers::sumWellRates(groupRein, schedule_, wellState, reportStepIdx, gasPos, /*isInjector*/false);
double gasProductionRate = WellGroupHelpers::sumWellSurfaceRates(groupRein, schedule_, wellState, reportStepIdx, gasPos, /*isInjector*/false);
double solventProductionRate = WellGroupHelpers::sumSolventRates(groupRein, schedule_, wellState, reportStepIdx, /*isInjector*/false);
solventProductionRate = comm_.sum(solventProductionRate);
@@ -1173,7 +1172,7 @@ getGuideRateValues(const Well& well) const
auto grval = data::GuideRateValue{};
const auto& wname = well.name();
if (!this->wellState().hasWellRates(wname)) {
if (!this->wellState().has(wname)) {
// No flow rates for 'wname' -- might be before well comes
// online (e.g., for the initial condition before simulation
// starts).
@@ -1333,7 +1332,7 @@ calculateAllGroupGuiderates(const int reportStepIdx) const
// group tree (FIELD group).
for (const auto& wname : schedule_.wellNames(reportStepIdx)) {
if (! (this->wellState().hasWellRates(wname) &&
if (! (this->wellState().has(wname) &&
this->guideRate_.has(wname)))
{
continue;
@@ -1528,7 +1527,7 @@ updateAndCommunicateGroupData(const int reportStepIdx,
// Set ALQ for off-process wells to zero
for (const auto& wname : schedule().wellNames(reportStepIdx)) {
const bool is_producer = schedule().getWell(wname, reportStepIdx).isProducer();
const bool not_on_this_process = well_state.wellMap().count(wname) == 0;
const bool not_on_this_process = !well_state.has(wname);
if (is_producer && not_on_this_process) {
well_state.setALQ(wname, 0.0);
}
opm/simulators/wells/BlackoilWellModel_impl.hpp
+2 -1
View File
@@ -1707,7 +1707,8 @@ namespace Opm {
auto& well_info = *local_parallel_well_info_[wellID];
const int num_perf_this_well = well_info.communication().sum(well_perf_data_[wellID].size());
auto& perf_data = this->wellState().perfData(wellID);
auto& ws = this->wellState().well(wellID);
auto& perf_data = ws.perf_data;
auto& perf_phase_rate = perf_data.phase_rates;
for (int perf = 0; perf < num_perf_this_well; ++perf) {
opm/simulators/wells/GasLiftGroupInfo.cpp
+2 -1
View File
@@ -300,7 +300,8 @@ GasLiftGroupInfo::
getProducerWellRates_(int well_index)
{
const auto& pu = this->phase_usage_;
const auto& wrate = this->well_state_.wellRates(well_index);
const auto& ws= this->well_state_.well(well_index);
const auto& wrate = ws.surface_rates;
const auto oil_rate = pu.phase_used[Oil]
? -wrate[pu.phase_pos[Oil]]
opm/simulators/wells/GasLiftSingleWellGeneric.cpp
+1 -1
View File
@@ -1235,7 +1235,7 @@ bool
GasLiftSingleWellGeneric::OptimizeState::
checkThpControl()
{
const int well_index = this->parent.well_state_.wellIndex(this->parent.well_name_);
const int well_index = this->parent.well_state_.index(this->parent.well_name_).value();
const Well::ProducerCMode& control_mode = this->parent.well_state_.well(well_index).production_cmode;
return control_mode == Well::ProducerCMode::THP;
}
opm/simulators/wells/GasLiftStage2.cpp
+3 -4
View File
@@ -404,11 +404,10 @@ GasLiftStage2::
getStdWellRates_(const WellInterfaceGeneric &well)
{
const int well_index = well.indexOfWell();
const auto& ws = this->well_state_.well(well_index);
const auto& pu = well.phaseUsage();
auto oil_rate =
-this->well_state_.wellRates(well_index)[pu.phase_pos[Oil]];
auto gas_rate =
-this->well_state_.wellRates(well_index)[pu.phase_pos[Gas]];
auto oil_rate = -ws.surface_rates[pu.phase_pos[Oil]];
auto gas_rate = -ws.surface_rates[pu.phase_pos[Gas]];
return {oil_rate, gas_rate};
}
opm/simulators/wells/MultisegmentWellEval.cpp
+7 -9
View File
@@ -1473,11 +1473,11 @@ updateThp(WellState& well_state,
static constexpr int Gas = BlackoilPhases::Vapour;
static constexpr int Oil = BlackoilPhases::Liquid;
static constexpr int Water = BlackoilPhases::Aqua;
auto& ws = well_state.well(baseif_.indexOfWell());
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
auto& well = well_state.well(baseif_.indexOfWell());
well.thp = 0;
ws.thp = 0;
return;
}
@@ -1486,18 +1486,16 @@ updateThp(WellState& well_state,
const PhaseUsage& pu = baseif_.phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Water ] ];
rates[ Water ] = ws.surface_rates[pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Oil ] ];
rates[ Oil ] = ws.surface_rates[pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Gas ] ];
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
auto& well = well_state.well(baseif_.indexOfWell());
const double bhp = well.bhp;
well.thp = this->calculateThpFromBhp(rates, bhp, rho, deferred_logger);
ws.thp = this->calculateThpFromBhp(rates, ws.bhp, rho, deferred_logger);
}
template<typename FluidSystem, typename Indices, typename Scalar>
@@ -1657,7 +1655,7 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
const double phase_rate = g_total * fractions[p];
segment_rates[seg*baseif_.numPhases() + p] = phase_rate;
if (seg == 0) { // top segment
well_state.wellRates(baseif_.indexOfWell())[p] = phase_rate;
ws.surface_rates[p] = phase_rate;
}
}
opm/simulators/wells/MultisegmentWellGeneric.cpp
+2 -2
View File
@@ -105,7 +105,7 @@ scaleSegmentRatesWithWellRates(WellState& well_state) const
auto& segment_rates = segments.rates;
for (int phase = 0; phase < baseif_.numPhases(); ++phase) {
const double unscaled_top_seg_rate = segment_rates[phase];
const double well_phase_rate = well_state.wellRates(baseif_.indexOfWell())[phase];
const double well_phase_rate = ws.surface_rates[phase];
if (std::abs(unscaled_top_seg_rate) > 1e-12)
{
for (int seg = 0; seg < numberOfSegments(); ++seg) {
@@ -120,7 +120,7 @@ scaleSegmentRatesWithWellRates(WellState& well_state) const
}
std::vector<double> perforation_rates(baseif_.numPhases() * baseif_.numPerfs(),0.0);
const double perf_phaserate_scaled = well_state.wellRates(baseif_.indexOfWell())[phase] / sumTw;
const double perf_phaserate_scaled = ws.surface_rates[phase] / sumTw;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
perforation_rates[baseif_.numPhases()* perf + phase] = baseif_.wellIndex()[perf] * perf_phaserate_scaled;
}
opm/simulators/wells/MultisegmentWell_impl.hpp
+5 -5
View File
@@ -279,14 +279,14 @@ namespace Opm
double total_rate = 0.0;
for (int phase = 0; phase < np; ++phase){
total_rate += well_state.wellRates(index_of_well_)[phase];
total_rate += ws.surface_rates[phase];
}
// for pressure controlled wells the well rates are the potentials
// if the rates are trivial we are most probably looking at the newly
// opened well and we therefore make the affort of computing the potentials anyway.
if (std::abs(total_rate) > 0) {
for (int phase = 0; phase < np; ++phase){
well_potentials[phase] = well_state.wellRates(index_of_well_)[phase];
well_potentials[phase] = ws.surface_rates[phase];
}
return;
}
@@ -367,7 +367,7 @@ namespace Opm
const int np = number_of_phases_;
const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
for (int phase = 0; phase < np; ++phase){
well_state_copy.wellRates(well_copy.index_of_well_)[phase] = sign * ws.well_potentials[phase];
ws.surface_rates[phase] = sign * ws.well_potentials[phase];
}
well_copy.scaleSegmentRatesWithWellRates(well_state_copy);
@@ -598,7 +598,7 @@ namespace Opm
};
auto& ws = well_state.well(this->index_of_well_);
auto& perf_data = well_state.perfData(this->index_of_well_);
auto& perf_data = ws.perf_data;
auto* connPI = perf_data.prod_index.data();
auto* wellPI = ws.productivity_index.data();
@@ -1267,7 +1267,7 @@ namespace Opm
// calculating the perforation rate for each perforation that belongs to this segment
const EvalWell seg_pressure = this->getSegmentPressure(seg);
auto& perf_data = well_state.perfData(this->index_of_well_);
auto& perf_data = ws.perf_data;
auto& perf_rates = perf_data.phase_rates;
auto& perf_press_state = perf_data.pressure;
for (const int perf : this->segment_perforations_[seg]) {
opm/simulators/wells/PerfData.cpp
+8 -5
View File
@@ -24,16 +24,15 @@ namespace Opm
{
PerfData::PerfData(std::size_t num_perf, bool injector_, const PhaseUsage& pu_arg)
: pu(pu_arg)
, injector(injector_)
PerfData::PerfData(std::size_t num_perf, bool injector_, std::size_t num_phases)
: injector(injector_)
, pressure(num_perf)
, rates(num_perf)
, phase_rates(num_perf * pu.num_phases)
, phase_rates(num_perf * num_phases)
, solvent_rates(num_perf)
, polymer_rates(num_perf)
, brine_rates(num_perf)
, prod_index(num_perf * pu.num_phases)
, prod_index(num_perf * num_phases)
, cell_index(num_perf)
, connection_transmissibility_factor(num_perf)
, satnum_id(num_perf)
@@ -50,6 +49,10 @@ std::size_t PerfData::size() const {
return this->pressure.size();
}
bool PerfData::empty() const {
return this->pressure.empty();
}
bool PerfData::try_assign(const PerfData& other) {
if (this->size() != other.size())
return false;
opm/simulators/wells/PerfData.hpp
+2 -4
View File
@@ -23,20 +23,18 @@
#include <vector>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm
{
class PerfData
{
private:
PhaseUsage pu;
bool injector;
public:
PerfData(std::size_t num_perf, bool injector_, const PhaseUsage& pu);
PerfData(std::size_t num_perf, bool injector_, std::size_t num_phases);
std::size_t size() const;
bool empty() const;
bool try_assign(const PerfData& other);
opm/simulators/wells/SingleWellState.cpp
+7 -1
View File
@@ -21,11 +21,14 @@
namespace Opm {
SingleWellState::SingleWellState(bool is_producer, std::size_t num_phases, double temp)
SingleWellState::SingleWellState(bool is_producer, std::size_t num_perf, std::size_t num_phases, double temp)
: producer(is_producer)
, temperature(temp)
, well_potentials(num_phases)
, productivity_index(num_phases)
, surface_rates(num_phases)
, reservoir_rates(num_phases)
, perf_data(num_perf, !is_producer, num_phases)
{}
@@ -49,6 +52,9 @@ void SingleWellState::shut() {
this->bhp = 0;
this->thp = 0;
this->status = Well::Status::SHUT;
std::fill(this->surface_rates.begin(), this->surface_rates.end(), 0);
std::fill(this->reservoir_rates.begin(), this->reservoir_rates.end(), 0);
std::fill(this->productivity_index.begin(), this->productivity_index.end(), 0);
}
void SingleWellState::stop() {
opm/simulators/wells/SingleWellState.hpp
+5 -1
View File
@@ -25,12 +25,13 @@
#include <opm/simulators/wells/SegmentState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
#include <opm/simulators/wells/PerfData.hpp>
namespace Opm {
class SingleWellState {
public:
SingleWellState(bool is_producer, std::size_t num_phases, double temp);
SingleWellState(bool is_producer, std::size_t num_perf, std::size_t num_phases, double temp);
Well::Status status{Well::Status::OPEN};
bool producer;
@@ -41,6 +42,9 @@ public:
double vaporized_oil_rate{0};
std::vector<double> well_potentials;
std::vector<double> productivity_index;
std::vector<double> surface_rates;
std::vector<double> reservoir_rates;
PerfData perf_data;
SegmentState segments;
Events events;
Well::InjectorCMode injection_cmode{Well::InjectorCMode::CMODE_UNDEFINED};
opm/simulators/wells/StandardWellEval.cpp
+24 -25
View File
@@ -255,11 +255,11 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
const int well_index = baseif_.indexOfWell();
const int np = baseif_.numPhases();
const auto& pu = baseif_.phaseUsage();
const auto& well = well_state.well(well_index);
const auto& ws = well_state.well(well_index);
// the weighted total well rate
double total_well_rate = 0.0;
for (int p = 0; p < np; ++p) {
total_well_rate += baseif_.scalingFactor(p) * well_state.wellRates(well_index)[p];
total_well_rate += baseif_.scalingFactor(p) * ws.surface_rates[p];
}
// Not: for the moment, the first primary variable for the injectors is not G_total. The injection rate
@@ -267,13 +267,13 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
if (baseif_.isInjector()) {
switch (baseif_.wellEcl().injectorType()) {
case InjectorType::WATER:
primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Water]];
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Water]];
break;
case InjectorType::GAS:
primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Gas]];
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Gas]];
break;
case InjectorType::OIL:
primary_variables_[WQTotal] = well_state.wellRates(well_index)[pu.phase_pos[Oil]];
primary_variables_[WQTotal] = ws.surface_rates[pu.phase_pos[Oil]];
break;
case InjectorType::MULTI:
// Not supported.
@@ -287,10 +287,10 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
if (std::abs(total_well_rate) > 0.) {
if constexpr (has_wfrac_variable) {
primary_variables_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * well_state.wellRates(well_index)[pu.phase_pos[Water]] / total_well_rate;
primary_variables_[WFrac] = baseif_.scalingFactor(pu.phase_pos[Water]) * ws.surface_rates[pu.phase_pos[Water]] / total_well_rate;
}
if constexpr (has_gfrac_variable) {
primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (well_state.wellRates(well_index)[pu.phase_pos[Gas]]
primary_variables_[GFrac] = baseif_.scalingFactor(pu.phase_pos[Gas]) * (ws.surface_rates[pu.phase_pos[Gas]]
- (Indices::enableSolvent ? well_state.solventWellRate(well_index) : 0.0) ) / total_well_rate ;
}
if constexpr (Indices::enableSolvent) {
@@ -343,7 +343,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
// BHP
primary_variables_[Bhp] = well.bhp;
primary_variables_[Bhp] = ws.bhp;
}
template<class FluidSystem, class Indices, class Scalar>
@@ -557,11 +557,11 @@ updateThp(WellState& well_state,
static constexpr int Gas = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Gas;
static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
auto& well = well_state.well(baseif_.indexOfWell());
auto& ws = well_state.well(baseif_.indexOfWell());
// When there is no vaild VFP table provided, we set the thp to be zero.
if (!baseif_.isVFPActive(deferred_logger) || baseif_.wellIsStopped()) {
well.thp = 0;
ws.thp = 0;
return;
}
@@ -570,20 +570,18 @@ updateThp(WellState& well_state,
const PhaseUsage& pu = baseif_.phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
rates[ Water ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Water ] ];
rates[ Water ] = ws.surface_rates[pu.phase_pos[ Water ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
rates[ Oil ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Oil ] ];
rates[ Oil ] = ws.surface_rates[pu.phase_pos[ Oil ] ];
}
if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
rates[ Gas ] = well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[ Gas ] ];
rates[ Gas ] = ws.surface_rates[pu.phase_pos[ Gas ] ];
}
const double bhp = well.bhp;
well.thp = this->calculateThpFromBhp(well_state,
ws.thp = this->calculateThpFromBhp(well_state,
rates,
bhp,
ws.bhp,
deferred_logger);
}
@@ -656,29 +654,29 @@ updateWellStateFromPrimaryVariables(WellState& well_state,
F[pu.phase_pos[Gas]] += F_solvent;
}
auto& well = well_state.well(baseif_.indexOfWell());
well.bhp = primary_variables_[Bhp];
auto& ws = well_state.well(baseif_.indexOfWell());
ws.bhp = primary_variables_[Bhp];
// calculate the phase rates based on the primary variables
// for producers, this is not a problem, while not sure for injectors here
if (baseif_.isProducer()) {
const double g_total = primary_variables_[WQTotal];
for (int p = 0; p < baseif_.numPhases(); ++p) {
well_state.wellRates(baseif_.indexOfWell())[p] = g_total * F[p];
ws.surface_rates[p] = g_total * F[p];
}
} else { // injectors
for (int p = 0; p < baseif_.numPhases(); ++p) {
well_state.wellRates(baseif_.indexOfWell())[p] = 0.0;
ws.surface_rates[p] = 0.0;
}
switch (baseif_.wellEcl().injectorType()) {
case InjectorType::WATER:
well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Water]] = primary_variables_[WQTotal];
ws.surface_rates[pu.phase_pos[Water]] = primary_variables_[WQTotal];
break;
case InjectorType::GAS:
well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Gas]] = primary_variables_[WQTotal];
ws.surface_rates[pu.phase_pos[Gas]] = primary_variables_[WQTotal];
break;
case InjectorType::OIL:
well_state.wellRates(baseif_.indexOfWell())[pu.phase_pos[Oil]] = primary_variables_[WQTotal];
ws.surface_rates[pu.phase_pos[Oil]] = primary_variables_[WQTotal];
break;
case InjectorType::MULTI:
// Not supported.
@@ -697,7 +695,8 @@ StandardWellEval<FluidSystem,Indices,Scalar>::
updateWellStateFromPrimaryVariablesPolyMW(WellState& well_state) const
{
if (baseif_.isInjector()) {
auto& perf_data = well_state.perfData(baseif_.indexOfWell());
auto& ws = well_state.well(baseif_.indexOfWell());
auto& perf_data = ws.perf_data;
auto& perf_water_velocity = perf_data.water_velocity;
auto& perf_skin_pressure = perf_data.skin_pressure;
for (int perf = 0; perf < baseif_.numPerfs(); ++perf) {
opm/simulators/wells/StandardWell_impl.hpp
+20 -15
View File
@@ -431,7 +431,7 @@ namespace Opm
const int np = number_of_phases_;
std::vector<RateVector> connectionRates = connectionRates_; // Copy to get right size.
auto& perf_data = well_state.perfData(this->index_of_well_);
auto& perf_data = ws.perf_data;
auto& perf_rates = perf_data.phase_rates;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
// Calculate perforation quantities.
@@ -546,8 +546,8 @@ namespace Opm
computePerfRateEval(intQuants, mob, bhp, Tw, perf, allow_cf,
cq_s, perf_dis_gas_rate, perf_vap_oil_rate, deferred_logger);
auto& perf_data = well_state.perfData(this->index_of_well_);
auto& ws = well_state.well(this->index_of_well_);
auto& perf_data = ws.perf_data;
if constexpr (has_polymer && Base::has_polymermw) {
if (this->isInjector()) {
// Store the original water flux computed from the reservoir quantities.
@@ -1225,7 +1225,7 @@ namespace Opm
}
// Compute the average pressure in each well block
const auto& perf_press = well_state.perfData(w).pressure;
const auto& perf_press = ws.perf_data.pressure;
auto p_above = this->parallel_well_info_.communicateAboveValues(ws.bhp,
perf_press.data(),
nperf);
@@ -1252,10 +1252,10 @@ namespace Opm
const int gaspos = gasCompIdx + perf * num_components_;
if (oilPresent) {
const double oilrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]); //in order to handle negative rates in producers
rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
if (oilrate > 0) {
const double gasrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
double rv = 0.0;
if (gasrate > 0) {
rv = oilrate / gasrate;
@@ -1278,9 +1278,9 @@ namespace Opm
const int oilpos = oilCompIdx + perf * num_components_;
if (gasPresent) {
rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
const double gasrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
const double gasrate = std::abs(ws.surface_rates[pu.phase_pos[Gas]]) - (has_solvent ? well_state.solventWellRate(w) : 0.0);
if (gasrate > 0) {
const double oilrate = std::abs(well_state.wellRates(w)[pu.phase_pos[Oil]]);
const double oilrate = std::abs(ws.surface_rates[pu.phase_pos[Oil]]);
double rs = 0.0;
if (oilrate > 0) {
rs = gasrate / oilrate;
@@ -1375,7 +1375,7 @@ namespace Opm
};
auto& ws = well_state.well(this->index_of_well_);
auto& perf_data = well_state.perfData(this->index_of_well_);
auto& perf_data = ws.perf_data;
auto* wellPI = ws.productivity_index.data();
auto* connPI = perf_data.prod_index.data();
@@ -1438,7 +1438,8 @@ namespace Opm
const int nperf = number_of_perforations_;
const int np = number_of_phases_;
std::vector<double> perfRates(b_perf.size(),0.0);
const auto& perf_data = well_state.perfData(this->index_of_well_);
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_data = ws.perf_data;
const auto& perf_rates_state = perf_data.phase_rates;
for (int perf = 0; perf < nperf; ++perf) {
@@ -1862,14 +1863,14 @@ namespace Opm
double total_rate = 0.0;
for (int phase = 0; phase < np; ++phase){
total_rate += well_state.wellRates(index_of_well_)[phase];
total_rate += ws.surface_rates[phase];
}
// for pressure controlled wells the well rates are the potentials
// if the rates are trivial we are most probably looking at the newly
// opened well and we therefore make the affort of computing the potentials anyway.
if (std::abs(total_rate) > 0) {
for (int phase = 0; phase < np; ++phase){
well_potentials[phase] = well_state.wellRates(index_of_well_)[phase];
well_potentials[phase] = ws.surface_rates[phase];
}
return;
}
@@ -1908,7 +1909,8 @@ namespace Opm
// other primary variables related to polymer injection
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
const auto& perf_data = well_state.perfData(this->index_of_well_);
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_data = ws.perf_data;
const auto& water_velocity = perf_data.water_velocity;
const auto& skin_pressure = perf_data.skin_pressure;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
@@ -2134,7 +2136,8 @@ namespace Opm
{
if constexpr (Base::has_polymermw) {
if (this->isInjector()) {
auto& perf_water_throughput = well_state.perfData(this->index_of_well_).water_throughput;
auto& ws = well_state.well(this->index_of_well_);
auto& perf_water_throughput = ws.perf_data.water_throughput;
for (int perf = 0; perf < number_of_perforations_; ++perf) {
const double perf_water_vel = this->primary_variables_[Bhp + 1 + perf];
// we do not consider the formation damage due to water flowing from reservoir into wellbore
@@ -2195,7 +2198,8 @@ namespace Opm
const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
this->resWell_[0][wat_vel_index] = eq_wat_vel.value();
const auto& perf_data = well_state.perfData(this->index_of_well_);
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_data = ws.perf_data;
const auto& perf_water_throughput = perf_data.water_throughput;
const double throughput = perf_water_throughput[perf];
const int pskin_index = Bhp + 1 + number_of_perforations_ + perf;
@@ -2257,7 +2261,8 @@ namespace Opm
const int wat_vel_index = Bhp + 1 + perf;
const EvalWell water_velocity = this->primary_variables_evaluation_[wat_vel_index];
if (water_velocity > 0.) { // injecting
const auto& perf_water_throughput = well_state.perfData(this->index_of_well_).water_throughput;
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_water_throughput = ws.perf_data.water_throughput;
const double throughput = perf_water_throughput[perf];
const EvalWell molecular_weight = wpolymermw(throughput, water_velocity, deferred_logger);
cq_s_polymw *= molecular_weight;
opm/simulators/wells/WellGroupHelpers.cpp
+100 -109
View File
@@ -58,6 +58,59 @@ namespace {
return {oilRate, gasRate, waterRate};
}
double sumWellPhaseRates(bool res_rates,
const Opm::Group& group,
const Opm::Schedule& schedule,
const Opm::WellState& wellState,
const int reportStepIdx,
const int phasePos,
const bool injector)
{
double rate = 0.0;
for (const std::string& groupName : group.groups()) {
const auto& groupTmp = schedule.getGroup(groupName, reportStepIdx);
const auto& gefac = groupTmp.getGroupEfficiencyFactor();
rate += gefac * sumWellPhaseRates(res_rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector);
}
for (const std::string& wellName : group.wells()) {
const auto& well_index = wellState.index(wellName);
if (!well_index.has_value())
continue;
if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once
{
continue;
}
const auto& wellEcl = schedule.getWell(wellName, reportStepIdx);
// only count producers or injectors
if ((wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector))
continue;
if (wellEcl.getStatus() == Opm::Well::Status::SHUT)
continue;
double factor = wellEcl.getEfficiencyFactor();
const auto& ws = wellState.well(well_index.value());
if (res_rates) {
const auto& well_rates = ws.reservoir_rates;
if (injector)
rate += factor * well_rates[phasePos];
else
rate -= factor * well_rates[phasePos];
} else {
const auto& well_rates = ws.surface_rates;
if (injector)
rate += factor * well_rates[phasePos];
else
rate -= factor * well_rates[phasePos];
}
}
return rate;
}
} // namespace Anonymous
namespace Opm
@@ -126,61 +179,15 @@ namespace WellGroupHelpers
schedule.getGroup(group.parent(), reportStepIdx), schedule, reportStepIdx, factor);
}
double sumWellPhaseRates(const WellContainer<std::vector<double>>& rates,
const Group& group,
const Schedule& schedule,
const WellState& wellState,
const int reportStepIdx,
const int phasePos,
const bool injector)
double sumWellSurfaceRates(const Group& group,
const Schedule& schedule,
const WellState& wellState,
const int reportStepIdx,
const int phasePos,
const bool injector)
{
double rate = 0.0;
for (const std::string& groupName : group.groups()) {
const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
const auto& gefac = groupTmp.getGroupEfficiencyFactor();
rate += gefac * sumWellPhaseRates(rates, groupTmp, schedule, wellState, reportStepIdx, phasePos, injector);
}
const auto& end = wellState.wellMap().end();
for (const std::string& wellName : group.wells()) {
const auto& it = wellState.wellMap().find(wellName);
if (it == end) // the well is not found
continue;
int well_index = it->second[0];
if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
{
continue;
}
const auto& wellEcl = schedule.getWell(wellName, reportStepIdx);
// only count producers or injectors
if ((wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector))
continue;
if (wellEcl.getStatus() == Well::Status::SHUT)
continue;
double factor = wellEcl.getEfficiencyFactor();
const auto& well_rates = rates[well_index];
if (injector)
rate += factor * well_rates[phasePos];
else
rate -= factor * well_rates[phasePos];
}
return rate;
}
double sumWellRates(const Group& group,
const Schedule& schedule,
const WellState& wellState,
const int reportStepIdx,
const int phasePos,
const bool injector)
{
return sumWellPhaseRates(wellState.wellRates(), group, schedule, wellState, reportStepIdx, phasePos, injector);
return sumWellPhaseRates(false, group, schedule, wellState, reportStepIdx, phasePos, injector);
}
double sumWellResRates(const Group& group,
@@ -190,8 +197,7 @@ namespace WellGroupHelpers
const int phasePos,
const bool injector)
{
return sumWellPhaseRates(
wellState.wellReservoirRates(), group, schedule, wellState, reportStepIdx, phasePos, injector);
return sumWellPhaseRates(true, group, schedule, wellState, reportStepIdx, phasePos, injector);
}
double sumSolventRates(const Group& group,
@@ -207,15 +213,13 @@ namespace WellGroupHelpers
const auto& gefac = groupTmp.getGroupEfficiencyFactor();
rate += gefac * sumSolventRates(groupTmp, schedule, wellState, reportStepIdx, injector);
}
const auto& end = wellState.wellMap().end();
for (const std::string& wellName : group.wells()) {
const auto& it = wellState.wellMap().find(wellName);
if (it == end) // the well is not found
const auto& well_index = wellState.index(wellName);
if (!well_index.has_value())
continue;
int well_index = it->second[0];
if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once
{
continue;
}
@@ -230,9 +234,9 @@ namespace WellGroupHelpers
double factor = wellEcl.getEfficiencyFactor();
if (injector)
rate += factor * wellState.solventWellRate(well_index);
rate += factor * wellState.solventWellRate(well_index.value());
else
rate -= factor * wellState.solventWellRate(well_index);
rate -= factor * wellState.solventWellRate(well_index.value());
}
return rate;
}
@@ -341,7 +345,7 @@ namespace WellGroupHelpers
if (individual_control || num_group_controlled_wells == 0) {
for (int phase = 0; phase < np; phase++) {
groupTargetReduction[phase]
+= subGroupEfficiency * sumWellRates(subGroup, schedule, wellStateNupcol, reportStepIdx, phase, isInjector);
+= subGroupEfficiency * sumWellSurfaceRates(subGroup, schedule, wellStateNupcol, reportStepIdx, phase, isInjector);
}
} else {
// The subgroup may participate in group control.
@@ -366,7 +370,7 @@ namespace WellGroupHelpers
if (individual_control || num_group_controlled_wells == 0) {
for (int phase = 0; phase < np; phase++) {
groupTargetReduction[phase]
+= subGroupEfficiency * sumWellRates(subGroup, schedule, wellStateNupcol, reportStepIdx, phase, isInjector);
+= subGroupEfficiency * sumWellSurfaceRates(subGroup, schedule, wellStateNupcol, reportStepIdx, phase, isInjector);
}
} else {
// The subgroup may participate in group control.
@@ -392,30 +396,28 @@ namespace WellGroupHelpers
if (wellTmp.getStatus() == Well::Status::SHUT)
continue;
const auto& end = wellState.wellMap().end();
const auto& it = wellState.wellMap().find(wellName);
if (it == end) // the well is not found
const auto& well_index = wellState.index(wellName);
if (!well_index.has_value())
continue;
int well_index = it->second[0];
if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once
{
continue;
}
const double efficiency = wellTmp.getEfficiencyFactor();
// add contributino from wells not under group control
const auto& ws = wellState.well(well_index);
const auto& ws_nupcol = wellStateNupcol.well(well_index.value());
const auto& ws = wellState.well(well_index.value());
if (isInjector) {
if (ws.injection_cmode != Well::InjectorCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
groupTargetReduction[phase] += wellStateNupcol.wellRates(well_index)[phase] * efficiency;
groupTargetReduction[phase] += ws_nupcol.surface_rates[phase] * efficiency;
}
} else {
if (ws.production_cmode != Well::ProducerCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
groupTargetReduction[phase] -= wellStateNupcol.wellRates(well_index)[phase] * efficiency;
groupTargetReduction[phase] -= ws_nupcol.surface_rates[phase] * efficiency;
}
}
}
@@ -466,29 +468,26 @@ namespace WellGroupHelpers
if (wellTmp.getStatus() == Well::Status::SHUT)
continue;
const auto& end = wellState.wellMap().end();
const auto& it = wellState.wellMap().find(wellName);
if (it == end) // the well is not found
const auto& well_index = wellState.index(wellName);
if (!well_index.has_value())
continue;
int well_index = it->second[0];
if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once
{
continue;
}
// scale rates
const auto& ws = wellState.well(well_index);
auto& ws = wellState.well(well_index.value());
if (isInjector) {
if (ws.injection_cmode == Well::InjectorCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
wellState.wellRates(well_index)[phase] *= scale;
ws.surface_rates[phase] *= scale;
}
} else {
if (ws.production_cmode == Well::ProducerCMode::GRUP)
for (int phase = 0; phase < np; phase++) {
wellState.wellRates(well_index)[phase] *= scale;
ws.surface_rates[phase] *= scale;
}
}
}
@@ -511,10 +510,10 @@ namespace WellGroupHelpers
const int np = wellState.numPhases();
double resv = 0.0;
for (int phase = 0; phase < np; ++phase) {
resv += sumWellPhaseRates(wellStateNupcol.wellReservoirRates(),
resv += sumWellPhaseRates(true,
group,
schedule,
wellState,
wellStateNupcol,
reportStepIdx,
phase,
/*isInjector*/ false);
@@ -536,10 +535,10 @@ namespace WellGroupHelpers
const int np = wellState.numPhases();
std::vector<double> resv(np, 0.0);
for (int phase = 0; phase < np; ++phase) {
resv[phase] = sumWellPhaseRates(wellStateNupcol.wellReservoirRates(),
resv[phase] = sumWellPhaseRates(true,
group,
schedule,
wellState,
wellStateNupcol,
reportStepIdx,
phase,
/*isInjector*/ true);
@@ -558,20 +557,19 @@ namespace WellGroupHelpers
updateWellRates(groupTmp, schedule, reportStepIdx, wellStateNupcol, wellState);
}
const int np = wellState.numPhases();
const auto& end = wellState.wellMap().end();
for (const std::string& wellName : group.wells()) {
std::vector<double> rates(np, 0.0);
const auto& it = wellState.wellMap().find(wellName);
if (it != end) { // the well is found on this node
int well_index = it->second[0];
const auto& well_index = wellState.index(wellName);
if (well_index.has_value()) { // the well is found on this node
const auto& wellTmp = schedule.getWell(wellName, reportStepIdx);
int sign = 1;
// production wellRates are negative. The users of currentWellRates uses the convention in
// opm-common that production and injection rates are positive.
if (!wellTmp.isInjector())
sign = -1;
const auto& ws = wellStateNupcol.well(well_index.value());
for (int phase = 0; phase < np; ++phase) {
rates[phase] = sign * wellStateNupcol.wellRates(well_index)[phase];
rates[phase] = sign * ws.surface_rates[phase];
}
}
wellState.setCurrentWellRates(wellName, rates);
@@ -592,8 +590,7 @@ namespace WellGroupHelpers
const int np = wellState.numPhases();
std::vector<double> rates(np, 0.0);
for (int phase = 0; phase < np; ++phase) {
rates[phase] = sumWellPhaseRates(
wellStateNupcol.wellRates(), group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false);
rates[phase] = sumWellPhaseRates(false, group, schedule, wellStateNupcol, reportStepIdx, phase, /*isInjector*/ false);
}
group_state.update_production_rates(group.name(), rates);
}
@@ -616,8 +613,7 @@ namespace WellGroupHelpers
std::vector<double> rein(np, 0.0);
for (int phase = 0; phase < np; ++phase) {
rein[phase] = sumWellPhaseRates(
wellStateNupcol.wellRates(), group, schedule, wellState, reportStepIdx, phase, /*isInjector*/ false);
rein[phase] = sumWellPhaseRates(false, group, schedule, wellStateNupcol, reportStepIdx, phase, /*isInjector*/ false);
}
// add import rate and substract consumption rate for group for gas
@@ -1371,19 +1367,16 @@ namespace WellGroupHelpers
if (wellTmp.getStatus() == Well::Status::SHUT)
continue;
const auto& end = wellState.wellMap().end();
const auto& it = wellState.wellMap().find(wellName);
if (it == end) // the well is not found
const auto& well_index = wellState.index(wellName);
if (!well_index.has_value()) // the well is not found
continue;
int well_index = it->second[0];
if (! wellState.wellIsOwned(well_index, wellName) ) // Only sum once
if (! wellState.wellIsOwned(well_index.value(), wellName) ) // Only sum once
{
continue;
}
const auto& ws = wellState.well(well_index);
const auto& ws = wellState.well(well_index.value());
// add contribution from wells unconditionally
for (int phase = 0; phase < np; phase++) {
pot[phase] += wefac * ws.well_potentials[phase];
@@ -1421,18 +1414,16 @@ namespace WellGroupHelpers
const Comm& comm,
GuideRate* guideRate)
{
const auto& end = wellState.wellMap().end();
for (const auto& well : schedule.getWells(reportStepIdx)) {
double oilpot = 0.0;
double gaspot = 0.0;
double waterpot = 0.0;
const auto& it = wellState.wellMap().find(well.name());
if (it != end && wellState.wellIsOwned(it->second[0], well.name()))
const auto& well_index = wellState.index(well.name());
if (well_index.has_value() && wellState.wellIsOwned(well_index.value(), well.name()))
{
// the well is found and owned
int well_index = it->second[0];
const auto& ws = wellState.well(well_index);
const auto& ws = wellState.well(well_index.value());
const auto& wpot = ws.well_potentials;
if (pu.phase_used[BlackoilPhases::Liquid] > 0)
oilpot = wpot[pu.phase_pos[BlackoilPhases::Liquid]];
opm/simulators/wells/WellGroupHelpers.hpp
+6 -14
View File
@@ -61,20 +61,12 @@ namespace WellGroupHelpers
const int reportStepIdx,
double& factor);
double sumWellPhaseRates(const WellContainer<std::vector<double>>& rates,
const Group& group,
const Schedule& schedule,
const WellState& wellState,
const int reportStepIdx,
const int phasePos,
const bool injector);
double sumWellRates(const Group& group,
const Schedule& schedule,
const WellState& wellState,
const int reportStepIdx,
const int phasePos,
const bool injector);
double sumWellSurfaceRates(const Group& group,
const Schedule& schedule,
const WellState& wellState,
const int reportStepIdx,
const int phasePos,
const bool injector);
double sumWellResRates(const Group& group,
const Schedule& schedule,
opm/simulators/wells/WellInterfaceFluidSystem.cpp
+30 -28
View File
@@ -67,17 +67,15 @@ calculateReservoirRates(WellState& well_state) const
const int fipreg = 0; // not considering the region for now
const int np = number_of_phases_;
auto& ws = well_state.well(this->index_of_well_);
std::vector<double> surface_rates(np, 0.0);
for (int p = 0; p < np; ++p) {
surface_rates[p] = well_state.wellRates(index_of_well_)[p];
surface_rates[p] = ws.surface_rates[p];
}
std::vector<double> voidage_rates(np, 0.0);
rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegionIdx_, surface_rates, voidage_rates);
for (int p = 0; p < np; ++p) {
well_state.wellReservoirRates(index_of_well_)[p] = voidage_rates[p];
}
ws.reservoir_rates = voidage_rates;
}
@@ -101,26 +99,26 @@ activeProductionConstraint(const WellState& well_state,
}
if (controls.hasControl(Well::ProducerCMode::ORAT) && currentControl != Well::ProducerCMode::ORAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Liquid]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
if (controls.oil_rate < current_rate)
return Well::ProducerCMode::ORAT;
}
if (controls.hasControl(Well::ProducerCMode::WRAT) && currentControl != Well::ProducerCMode::WRAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Aqua]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
if (controls.water_rate < current_rate)
return Well::ProducerCMode::WRAT;
}
if (controls.hasControl(Well::ProducerCMode::GRAT) && currentControl != Well::ProducerCMode::GRAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Vapour]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]];
if (controls.gas_rate < current_rate)
return Well::ProducerCMode::GRAT;
}
if (controls.hasControl(Well::ProducerCMode::LRAT) && currentControl != Well::ProducerCMode::LRAT) {
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Liquid]];
current_rate -= well_state.wellRates(well_index)[pu.phase_pos[BlackoilPhases::Aqua]];
double current_rate = -ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]];
current_rate -= ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]];
if (controls.liquid_rate < current_rate)
return Well::ProducerCMode::LRAT;
}
@@ -128,13 +126,13 @@ activeProductionConstraint(const WellState& well_state,
if (controls.hasControl(Well::ProducerCMode::RESV) && currentControl != Well::ProducerCMode::RESV) {
double current_rate = 0.0;
if (pu.phase_used[BlackoilPhases::Aqua])
current_rate -= well_state.wellReservoirRates(well_index)[pu.phase_pos[BlackoilPhases::Aqua]];
current_rate -= ws.reservoir_rates[pu.phase_pos[BlackoilPhases::Aqua]];
if (pu.phase_used[BlackoilPhases::Liquid])
current_rate -= well_state.wellReservoirRates(well_index)[pu.phase_pos[BlackoilPhases::Liquid]];
current_rate -= ws.reservoir_rates[pu.phase_pos[BlackoilPhases::Liquid]];
if (pu.phase_used[BlackoilPhases::Vapour])
current_rate -= well_state.wellReservoirRates(well_index)[pu.phase_pos[BlackoilPhases::Vapour]];
current_rate -= ws.reservoir_rates[pu.phase_pos[BlackoilPhases::Vapour]];
if (controls.prediction_mode && controls.resv_rate < current_rate)
return Well::ProducerCMode::RESV;
@@ -203,17 +201,17 @@ activeInjectionConstraint(const WellState& well_state,
switch (injectorType) {
case InjectorType::WATER:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
current_rate = ws.surface_rates[ pu.phase_pos[BlackoilPhases::Aqua] ];
break;
}
case InjectorType::OIL:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
current_rate = ws.surface_rates[ pu.phase_pos[BlackoilPhases::Liquid] ];
break;
}
case InjectorType::GAS:
{
current_rate = well_state.wellRates(well_index)[ pu.phase_pos[BlackoilPhases::Vapour] ];
current_rate = ws.surface_rates[ pu.phase_pos[BlackoilPhases::Vapour] ];
break;
}
default:
@@ -228,13 +226,13 @@ activeInjectionConstraint(const WellState& well_state,
{
double current_rate = 0.0;
if( pu.phase_used[BlackoilPhases::Aqua] )
current_rate += well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Aqua] ];
current_rate += ws.reservoir_rates[ pu.phase_pos[BlackoilPhases::Aqua] ];
if( pu.phase_used[BlackoilPhases::Liquid] )
current_rate += well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Liquid] ];
current_rate += ws.reservoir_rates[ pu.phase_pos[BlackoilPhases::Liquid] ];
if( pu.phase_used[BlackoilPhases::Vapour] )
current_rate += well_state.wellReservoirRates(well_index)[ pu.phase_pos[BlackoilPhases::Vapour] ];
current_rate += ws.reservoir_rates[ pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.reservoir_rate < current_rate)
return Well::InjectorCMode::RESV;
@@ -317,6 +315,7 @@ checkGroupConstraintsInj(const Group& group,
std::vector<double> resv_coeff(phaseUsage().num_phases, 1.0);
rateConverter_.calcInjCoeff(0, pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
const auto& ws = well_state.well(this->index_of_well_);
// Call check for the well's injection phase.
return WellGroupHelpers::checkGroupConstraintsInj(name(),
well_ecl_.groupName(),
@@ -325,7 +324,7 @@ checkGroupConstraintsInj(const Group& group,
group_state,
current_step_,
guide_rate_,
well_state.wellRates(index_of_well_).data(),
ws.surface_rates.data(),
injectionPhase,
phaseUsage(),
efficiencyFactor,
@@ -350,6 +349,7 @@ checkGroupConstraintsProd(const Group& group,
std::vector<double> resv_coeff(this->phaseUsage().num_phases, 1.0);
rateConverter_.calcCoeff(0, pvtRegionIdx_, resv_coeff); // FIPNUM region 0 here, should use FIPNUM from WELSPECS.
const auto& ws = well_state.well(this->index_of_well_);
return WellGroupHelpers::checkGroupConstraintsProd(name(),
well_ecl_.groupName(),
group,
@@ -357,7 +357,7 @@ checkGroupConstraintsProd(const Group& group,
group_state,
current_step_,
guide_rate_,
well_state.wellRates(index_of_well_).data(),
ws.surface_rates.data(),
phaseUsage(),
efficiencyFactor,
schedule,
@@ -400,7 +400,7 @@ checkGroupConstraints(WellState& well_state,
ws.injection_cmode = Well::InjectorCMode::GRUP;
const int np = well_state.numPhases();
for (int p = 0; p<np; ++p) {
well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
ws.surface_rates[p] *= group_constraint.second;
}
}
return group_constraint.first;
@@ -428,7 +428,7 @@ checkGroupConstraints(WellState& well_state,
ws.production_cmode = Well::ProducerCMode::GRUP;
const int np = well_state.numPhases();
for (int p = 0; p<np; ++p) {
well_state.wellRates(index_of_well_)[p] *= group_constraint.second;
ws.surface_rates[p] *= group_constraint.second;
}
}
return group_constraint.first;
@@ -698,7 +698,7 @@ updateWellTestStateEconomic(const WellState& well_state,
if (quantity_limit == WellEconProductionLimits::QuantityLimit::POTN)
rate_limit_violated = checkRateEconLimits(econ_production_limits, ws.well_potentials.data(), deferred_logger);
else {
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellRates(index_of_well_).data(), deferred_logger);
rate_limit_violated = checkRateEconLimits(econ_production_limits, ws.surface_rates.data(), deferred_logger);
}
}
@@ -854,7 +854,8 @@ checkMaxRatioLimitCompletions(const WellState& well_state,
double max_ratio_completion = 0;
const int np = number_of_phases_;
const auto& perf_data = well_state.perfData(this->index_of_well_);
const auto& ws = well_state.well(this->index_of_well_);
const auto& perf_data = ws.perf_data;
const auto& perf_phase_rates = perf_data.phase_rates;
// look for the worst_offending_completion
for (const auto& completion : completions_) {
@@ -899,9 +900,9 @@ checkMaxRatioLimitWell(const WellState& well_state,
const int np = number_of_phases_;
std::vector<double> well_rates(np, 0.0);
const auto& ws = well_state.well(this->index_of_well_);
for (int p = 0; p < np; ++p) {
well_rates[p] = well_state.wellRates(index_of_well_)[p];
well_rates[p] = ws.surface_rates[p];
}
const double well_ratio = ratioFunc(well_rates, phaseUsage());
@@ -1099,7 +1100,8 @@ getGroupProductionTargetRate(const Group& group,
}
// Avoid negative target rates coming from too large local reductions.
const double target_rate = std::max(0.0, target / efficiencyFactor);
const auto& rates = well_state.wellRates(index_of_well_);
const auto& ws = well_state.well(this->index_of_well_);
const auto& rates = ws.surface_rates;
const auto current_rate = -tcalc.calcModeRateFromRates(rates); // Switch sign since 'rates' are negative for producers.
double scale = 1.0;
if (current_rate > 1e-14)
opm/simulators/wells/WellInterface_impl.hpp
+39 -48
View File
@@ -610,7 +610,7 @@ namespace Opm
if (this->wellIsStopped()) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = 0.0;
ws.surface_rates[p] = 0;
}
ws.thp = 0;
return;
@@ -647,7 +647,7 @@ namespace Opm
switch(current) {
case Well::InjectorCMode::RATE:
{
well_state.wellRates(well_index)[phasePos] = controls.surface_rate;
ws.surface_rates[phasePos] = controls.surface_rate;
break;
}
@@ -656,16 +656,13 @@ namespace Opm
std::vector<double> convert_coeff(this->number_of_phases_, 1.0);
this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
const double coeff = convert_coeff[phasePos];
well_state.wellRates(well_index)[phasePos] = controls.reservoir_rate/coeff;
ws.surface_rates[phasePos] = controls.reservoir_rate/coeff;
break;
}
case Well::InjectorCMode::THP:
{
std::vector<double> rates(3, 0.0);
for (int p = 0; p<np; ++p) {
rates[p] = well_state.wellRates(well_index)[p];
}
auto rates = ws.surface_rates;
double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
ws.bhp = bhp;
@@ -673,11 +670,9 @@ namespace Opm
// we try to provide a better intial well rate
// using the well potentials
double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = ws.well_potentials[p];
}
}
if (total_rate <= 0.0)
ws.surface_rates = ws.well_potentials;
break;
}
case Well::InjectorCMode::BHP:
@@ -685,16 +680,14 @@ namespace Opm
ws.bhp = controls.bhp_limit;
double total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate += well_state.wellRates(well_index)[p];
total_rate += ws.surface_rates[p];
}
// if the total rates are negative or zero
// we try to provide a better intial well rate
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = ws.well_potentials[p];
}
}
if (total_rate <= 0.0)
ws.surface_rates = ws.well_potentials;
break;
}
case Well::InjectorCMode::GRUP:
@@ -712,7 +705,7 @@ namespace Opm
efficiencyFactor,
deferred_logger);
if (target)
well_state.wellRates(well_index)[phasePos] = *target;
ws.surface_rates[phasePos] = *target;
break;
}
case Well::InjectorCMode::CMODE_UNDEFINED:
@@ -730,19 +723,19 @@ namespace Opm
switch (current) {
case Well::ProducerCMode::ORAT:
{
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[Oil] ];
double current_rate = -ws.surface_rates[ pu.phase_pos[Oil] ];
// for trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.oil_rate/current_rate;
ws.surface_rates[p] *= controls.oil_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Oil]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.oil_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.oil_rate/control_fraction;
}
}
}
@@ -750,19 +743,19 @@ namespace Opm
}
case Well::ProducerCMode::WRAT:
{
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[Water] ];
double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ];
// for trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.water_rate/current_rate;
ws.surface_rates[p] *= controls.water_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Water]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.water_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.water_rate/control_fraction;
}
}
}
@@ -770,19 +763,19 @@ namespace Opm
}
case Well::ProducerCMode::GRAT:
{
double current_rate = -well_state.wellRates(well_index)[pu.phase_pos[Gas] ];
double current_rate = -ws.surface_rates[pu.phase_pos[Gas] ];
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.gas_rate/current_rate;
ws.surface_rates[p] *= controls.gas_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Gas]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.gas_rate/control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.gas_rate/control_fraction;
}
}
}
@@ -792,20 +785,20 @@ namespace Opm
}
case Well::ProducerCMode::LRAT:
{
double current_rate = -well_state.wellRates(well_index)[ pu.phase_pos[Water] ]
- well_state.wellRates(well_index)[ pu.phase_pos[Oil] ];
double current_rate = -ws.surface_rates[ pu.phase_pos[Water] ]
- ws.surface_rates[ pu.phase_pos[Oil] ];
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (current_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.liquid_rate/current_rate;
ws.surface_rates[p] *= controls.liquid_rate/current_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
double control_fraction = fractions[pu.phase_pos[Water]] + fractions[pu.phase_pos[Oil]];
if (control_fraction != 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.liquid_rate / control_fraction;
ws.surface_rates[p] = - fractions[p] * controls.liquid_rate / control_fraction;
}
}
}
@@ -821,19 +814,19 @@ namespace Opm
this->rateConverter_.calcCoeff(/*fipreg*/ 0, this->pvtRegionIdx_, convert_coeff);
double total_res_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_res_rate -= well_state.wellRates(well_index)[p] * convert_coeff[p];
total_res_rate -= ws.surface_rates[p] * convert_coeff[p];
}
if (controls.prediction_mode) {
// or trivial rates or opposite direction we don't just scale the rates
// but use either the potentials or the mobility ratio to initial the well rates
if (total_res_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= controls.resv_rate/total_res_rate;
ws.surface_rates[p] *= controls.resv_rate/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
ws.surface_rates[p] = - fractions[p] * controls.resv_rate / convert_coeff[p];
}
}
} else {
@@ -854,12 +847,12 @@ namespace Opm
// but use either the potentials or the mobility ratio to initial the well rates
if (total_res_rate > 0.0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= target/total_res_rate;
ws.surface_rates[p] *= target/total_res_rate;
}
} else {
const std::vector<double> fractions = initialWellRateFractions(ebos_simulator, well_state);
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = - fractions[p] * target / convert_coeff[p];
ws.surface_rates[p] = - fractions[p] * target / convert_coeff[p];
}
}
@@ -871,24 +864,21 @@ namespace Opm
ws.bhp = controls.bhp_limit;
double total_rate = 0.0;
for (int p = 0; p<np; ++p) {
total_rate -= well_state.wellRates(well_index)[p];
total_rate -= ws.surface_rates[p];
}
// if the total rates are negative or zero
// we try to provide a better intial well rate
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = -ws.well_potentials[p];
ws.surface_rates[p] = -ws.well_potentials[p];
}
}
break;
}
case Well::ProducerCMode::THP:
{
std::vector<double> rates(3, 0.0);
for (int p = 0; p<np; ++p) {
rates[p] = well_state.wellRates(well_index)[p];
}
auto rates = ws.surface_rates;
double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
ws.bhp = bhp;
@@ -898,7 +888,7 @@ namespace Opm
double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = -ws.well_potentials[p];
ws.surface_rates[p] = -ws.well_potentials[p];
}
}
break;
@@ -918,7 +908,7 @@ namespace Opm
// we don't want to scale with zero and get zero rates.
if (scale > 0) {
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] *= scale;
ws.surface_rates[p] *= scale;
}
}
break;
@@ -989,9 +979,10 @@ namespace Opm
{
// Check if the rates of this well only are single-phase, do nothing
// if more than one nonzero rate.
auto& ws = well_state.well(this->index_of_well_);
int nonzero_rate_index = -1;
for (int p = 0; p < this->number_of_phases_; ++p) {
if (well_state.wellRates(this->index_of_well_)[p] != 0.0) {
if (ws.surface_rates[p] != 0.0) {
if (nonzero_rate_index == -1) {
nonzero_rate_index = p;
} else {
@@ -1009,12 +1000,12 @@ namespace Opm
std::vector<double> well_q_s = computeCurrentWellRates(ebosSimulator, deferred_logger);
// 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 double initial_nonzero_rate = ws.surface_rates[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 = this->flowPhaseToEbosCompIdx(p);
double& rate = well_state.wellRates(this->index_of_well_)[p];
double& rate = ws.surface_rates[p];
rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]);
}
}
opm/simulators/wells/WellState.cpp
+59 -118
View File
@@ -39,31 +39,17 @@ void WellState::base_init(const std::vector<double>& cellPressures,
const SummaryState& summary_state)
{
// clear old name mapping
this->wellMap_.clear();
this->perfdata.clear();
this->parallel_well_info_.clear();
this->wellrates_.clear();
this->wells_.clear();
{
// const int nw = wells->number_of_wells;
const int nw = wells_ecl.size();
// const int np = wells->number_of_phases;
int connpos = 0;
for (int w = 0; w < nw; ++w) {
const Well& well = wells_ecl[w];
// Initialize bhp(), thp(), wellRates(), temperature().
initSingleWell(cellPressures, w, well, well_perf_data[w], parallel_well_info[w], summary_state);
// Setup wellname -> well index mapping.
const int num_perf_this_well = well_perf_data[w].size();
std::string name = well.name();
assert( !name.empty() );
mapentry_t& wellMapEntry = wellMap_[name];
wellMapEntry[ 0 ] = w;
wellMapEntry[ 1 ] = connpos;
wellMapEntry[ 2 ] = num_perf_this_well;
connpos += num_perf_this_well;
initSingleWell(cellPressures, well, well_perf_data[w], parallel_well_info[w], summary_state);
}
}
}
@@ -73,11 +59,10 @@ void WellState::base_init(const std::vector<double>& cellPressures,
void WellState::initSingleWell(const std::vector<double>& cellPressures,
const int w,
const Well& well,
const std::vector<PerforationData>& well_perf_data,
const ParallelWellInfo* well_info,
const SummaryState& summary_state)
const Well& well,
const std::vector<PerforationData>& well_perf_data,
const ParallelWellInfo* well_info,
const SummaryState& summary_state)
{
assert(well.isInjector() || well.isProducer());
@@ -87,13 +72,10 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
const int np = pu.num_phases;
double temp = well.isInjector() ? well.injectionControls(summary_state).temperature : 273.15 + 15.56;
auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), static_cast<std::size_t>(np), temp});
this->parallel_well_info_.add(well.name(), well_info);
this->wellrates_.add(well.name(), std::vector<double>(np, 0));
const int num_perf_this_well = well_info->communication().sum(well_perf_data.size());
this->perfdata.add(well.name(), PerfData{static_cast<std::size_t>(num_perf_this_well), well.isInjector(), this->phase_usage_});
auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), well_perf_data.size(), static_cast<std::size_t>(np), temp});
if ( num_perf_this_well == 0 )
if ( ws.perf_data.empty())
return;
const auto inj_controls = well.isInjector() ? well.injectionControls(summary_state) : Well::InjectionControls(0);
@@ -142,7 +124,7 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
// (producer) or RATE (injector).
// Otherwise, we cannot set the correct
// value here and initialize to zero rate.
auto& rates = this->wellrates_[w];
auto& rates = ws.surface_rates;
if (well.isInjector()) {
if (inj_controls.cmode == Well::InjectorCMode::RATE) {
switch (inj_controls.injector_type) {
@@ -247,7 +229,6 @@ void WellState::init(const std::vector<double>& cellPressures,
nperf += wpd.size();
}
well_reservoir_rates_.clear();
{
const auto& wg_events = schedule[report_step].wellgroup_events();
for (const auto& ecl_well : wells_ecl) {
@@ -261,11 +242,10 @@ void WellState::init(const std::vector<double>& cellPressures,
// Initialize perfphaserates_ to well
// rates divided by the number of perforations.
const auto& ecl_well = wells_ecl[w];
const auto& wname = ecl_well.name();
const auto& well_info = this->wellMap().at(wname);
const int num_perf_this_well = well_info[2];
auto& ws = this->well(w);
auto& perf_data = ws.perf_data;
const int num_perf_this_well = perf_data.size();
const int global_num_perf_this_well = ecl_well.getConnections().num_open();
auto& perf_data = this->perfData(w);
const auto& perf_input = well_perf_data[w];
for (int perf = 0; perf < num_perf_this_well; ++perf) {
@@ -276,13 +256,11 @@ void WellState::init(const std::vector<double>& cellPressures,
if (wells_ecl[w].getStatus() == Well::Status::OPEN) {
for (int p = 0; p < this->numPhases(); ++p) {
perf_data.phase_rates[this->numPhases()*perf + p] = wellRates(w)[p] / double(global_num_perf_this_well);
perf_data.phase_rates[this->numPhases()*perf + p] = ws.surface_rates[p] / double(global_num_perf_this_well);
}
}
perf_data.pressure[perf] = cellPressures[well_perf_data[w][perf].cell_index];
}
this->well_reservoir_rates_.add(wname, std::vector<double>(np, 0));
}
for (int w = 0; w < nw; ++w) {
@@ -316,20 +294,16 @@ void WellState::init(const std::vector<double>& cellPressures,
// intialize wells that have been there before
// order may change so the mapping is based on the well name
if (prevState && !prevState->wellMap().empty()) {
auto end = prevState->wellMap().end();
if (prevState && prevState->size() > 0) {
for (int w = 0; w < nw; ++w) {
const Well& well = wells_ecl[w];
if (well.getStatus() == Well::Status::SHUT) {
continue;
}
auto& new_well = this->well(w);
auto it = prevState->wellMap().find(well.name());
if (it != end)
{
const int oldIndex = it->second[ 0 ];
const auto& prev_well = prevState->well(oldIndex);
const auto& old_index = prevState->index(well.name());
if (old_index.has_value()) {
const auto& prev_well = prevState->well(old_index.value());
new_well.init_timestep(prev_well);
@@ -348,21 +322,13 @@ void WellState::init(const std::vector<double>& cellPressures,
new_well.production_cmode = prev_well.production_cmode;
}
wellRates(w) = prevState->wellRates(oldIndex);
wellReservoirRates(w) = prevState->wellReservoirRates(oldIndex);
new_well.surface_rates = prev_well.surface_rates;
new_well.reservoir_rates = prev_well.reservoir_rates;
new_well.well_potentials = prev_well.well_potentials;
// perfPhaseRates
const int num_perf_old_well = (*it).second[ 2 ];
const auto new_iter = this->wellMap().find(well.name());
if (new_iter == this->wellMap().end()) {
throw std::logic_error {
well.name() + " is not in internal well map - "
"Bug in WellState"
};
}
const int num_perf_this_well = new_iter->second[2];
const int num_perf_old_well = prev_well.perf_data.size();
const int num_perf_this_well = new_well.perf_data.size();
const bool global_num_perf_same = (num_perf_this_well == num_perf_old_well);
// copy perforation rates when the number of
@@ -371,16 +337,16 @@ void WellState::init(const std::vector<double>& cellPressures,
// number of perforations.
if (global_num_perf_same)
{
auto& perf_data = this->perfData(w);
const auto& prev_perf_data = prevState->perfData(w);
auto& perf_data = new_well.perf_data;
const auto& prev_perf_data = prev_well.perf_data;
perf_data.try_assign( prev_perf_data );
} else {
const int global_num_perf_this_well = well.getConnections().num_open();
auto& perf_data = this->perfData(w);
auto& perf_data = new_well.perf_data;
auto& target_rates = perf_data.phase_rates;
for (int perf_index = 0; perf_index < num_perf_this_well; perf_index++) {
for (int p = 0; p < np; ++p) {
target_rates[perf_index*np + p] = wellRates(w)[p] / double(global_num_perf_this_well);
target_rates[perf_index*np + p] = new_well.surface_rates[p] / double(global_num_perf_this_well);
}
}
}
@@ -466,18 +432,19 @@ WellState::report(const int* globalCellIdxMap,
const auto& pu = this->phaseUsage();
data::Wells res;
for( const auto& [wname, winfo]: this->wellMap() ) {
const auto well_index = winfo[ 0 ];
for( std::size_t well_index = 0; well_index < this->size(); well_index++) {
const auto& ws = this->well(well_index);
if ((ws.status == Well::Status::SHUT) && !wasDynamicallyClosed(well_index))
{
continue;
}
const auto& reservoir_rates = this->well_reservoir_rates_[well_index];
const auto& reservoir_rates = ws.reservoir_rates;
const auto& well_potentials = ws.well_potentials;
const auto& wpi = ws.productivity_index;
const auto& wv = this->wellRates(well_index);
const auto& wv = ws.surface_rates;
const auto& wname = this->name(well_index);
data::Well well;
well.bhp = ws.bhp;
@@ -566,7 +533,7 @@ void WellState::reportConnections(std::vector<data::Connection>& connections,
const int* globalCellIdxMap) const
{
using rt = data::Rates::opt;
const auto& perf_data = this->perfData(well_index);
const auto& perf_data = this->well(well_index).perf_data;
const int num_perf_well = perf_data.size();
connections.resize(num_perf_well);
const auto& perf_rates = perf_data.rates;
@@ -579,8 +546,6 @@ void WellState::reportConnections(std::vector<data::Connection>& connections,
connection.reservoir_rate = perf_rates[i];
connection.trans_factor = perf_data.connection_transmissibility_factor[i];
}
assert(num_perf_well == int(connections.size()));
const int np = pu.num_phases;
size_t local_comp_index = 0;
@@ -623,7 +588,6 @@ void WellState::reportConnections(std::vector<data::Connection>& connections,
++local_comp_index;
}
assert(local_comp_index == this->perfdata[well_index].size());
}
void WellState::initWellStateMSWell(const std::vector<Well>& wells_ecl,
@@ -677,30 +641,24 @@ void WellState::initWellStateMSWell(const std::vector<Well>& wells_ecl,
}
auto& perf_data = this->perfData(w);
auto& perf_data = ws.perf_data;
// for the seg_rates_, now it becomes a recursive solution procedure.
{
// make sure the information from wells_ecl consistent with wells
assert((n_activeperf == this->wellMap().at(well_ecl.name())[2]) &&
"Inconsistent number of reservoir connections in well");
if (pu.phase_used[Gas]) {
auto& perf_rates = perf_data.phase_rates;
const int gaspos = pu.phase_pos[Gas];
// scale the phase rates for Gas to avoid too bad initial guess for gas fraction
// it will probably benefit the standard well too, while it needs to be justified
// TODO: to see if this strategy can benefit StandardWell too
// TODO: it might cause big problem for gas rate control or if there is a gas rate limit
// maybe the best way is to initialize the fractions first then get the rates
for (int perf = 0; perf < n_activeperf; perf++)
perf_rates[perf*np + gaspos] *= 100;
}
const auto& perf_rates = perf_data.phase_rates;
std::vector<double> perforation_rates(perf_rates.begin(), perf_rates.end());
calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates);
if (pu.phase_used[Gas]) {
auto& perf_rates = perf_data.phase_rates;
const int gaspos = pu.phase_pos[Gas];
// scale the phase rates for Gas to avoid too bad initial guess for gas fraction
// it will probably benefit the standard well too, while it needs to be justified
// TODO: to see if this strategy can benefit StandardWell too
// TODO: it might cause big problem for gas rate control or if there is a gas rate limit
// maybe the best way is to initialize the fractions first then get the rates
for (int perf = 0; perf < n_activeperf; perf++)
perf_rates[perf*np + gaspos] *= 100;
}
const auto& perf_rates = perf_data.phase_rates;
std::vector<double> perforation_rates(perf_rates.begin(), perf_rates.end());
calculateSegmentRates(segment_inlets, segment_perforations, perforation_rates, np, 0 /* top segment */, ws.segments.rates);
// for the segment pressure, the segment pressure is the same with the first perforation belongs to the segment
// if there is no perforation associated with this segment, it uses the pressure from the outlet segment
// which requres the ordering is successful
@@ -783,21 +741,24 @@ WellState::calculateSegmentRates(const std::vector<std::vector<int>>& segment_in
double WellState::solventWellRate(const int w) const
{
const auto& perf_data = this->perfData(w);
auto& ws = this->well(w);
const auto& perf_data = ws.perf_data;
const auto& perf_rates_solvent = perf_data.solvent_rates;
return parallel_well_info_[w]->sumPerfValues(perf_rates_solvent.begin(), perf_rates_solvent.end());
}
double WellState::polymerWellRate(const int w) const
{
const auto& perf_data = this->perfData(w);
auto& ws = this->well(w);
const auto& perf_data = ws.perf_data;
const auto& perf_rates_polymer = perf_data.polymer_rates;
return parallel_well_info_[w]->sumPerfValues(perf_rates_polymer.begin(), perf_rates_polymer.end());
}
double WellState::brineWellRate(const int w) const
{
const auto& perf_data = this->perfData(w);
auto& ws = this->well(w);
const auto& perf_data = ws.perf_data;
const auto& perf_rates_brine = perf_data.brine_rates;
return parallel_well_info_[w]->sumPerfValues(perf_rates_brine.begin(), perf_rates_brine.end());
}
@@ -813,18 +774,8 @@ void WellState::shutWell(int well_index)
{
auto& ws = this->well(well_index);
ws.shut();
const int np = numPhases();
this->wellrates_[well_index].assign(np, 0);
auto& resv = this->well_reservoir_rates_[well_index];
auto& wpi = ws.productivity_index;
for (int p = 0; p < np; ++p) {
resv[p] = 0.0;
wpi[p] = 0.0;
}
auto& perf_data = this->perfData(well_index);
auto& perf_data = ws.perf_data;
auto& connpi = perf_data.prod_index;
connpi.assign(connpi.size(), 0);
}
@@ -961,21 +912,11 @@ bool WellState::wellIsOwned(std::size_t well_index,
bool WellState::wellIsOwned(const std::string& wellName) const
{
const auto& it = this->wellMap_.find( wellName );
if (it == this->wellMap_.end()) {
const auto& well_index = this->index(wellName);
if (!well_index.has_value())
OPM_THROW(std::logic_error, "Could not find well " << wellName << " in well map");
}
const int well_index = it->second[0];
return wellIsOwned(well_index, wellName);
}
int WellState::wellIndex(const std::string& wellName) const
{
const auto& it = this->wellMap_.find( wellName );
if (it == this->wellMap_.end()) {
OPM_THROW(std::logic_error, "Could not find well " << wellName << " in well map");
}
return it->second[0];
return wellIsOwned(well_index.value(), wellName);
}
int WellState::numSegments(const int well_id) const
@@ -1006,16 +947,16 @@ void WellState::updateWellsDefaultALQ( const std::vector<Well>& wells_ecl )
void WellState::resetConnectionTransFactors(const int well_index,
const std::vector<PerforationData>& new_perf_data)
{
if (this->perfdata[well_index].size() != new_perf_data.size()) {
auto& ws = this->well(well_index);
auto& perf_data = ws.perf_data;
if (perf_data.size() != new_perf_data.size()) {
throw std::invalid_argument {
"Size mismatch for perforation data in well "
+ std::to_string(well_index)
};
}
auto& perf_data = this->perfData(well_index);
for (std::size_t conn_index = 0; conn_index < new_perf_data.size(); conn_index++) {
if (perf_data.cell_index[conn_index] != static_cast<std::size_t>(new_perf_data[conn_index].cell_index)) {
throw std::invalid_argument {
"Cell index mismatch in connection "
opm/simulators/wells/WellState.hpp
+7 -51
View File
@@ -52,9 +52,6 @@ class Schedule;
class WellState
{
public:
using mapentry_t = std::array<int, 3>;
using WellMapType = std::map<std::string, mapentry_t>;
static const uint64_t event_mask = ScheduleEvents::WELL_STATUS_CHANGE + ScheduleEvents::PRODUCTION_UPDATE + ScheduleEvents::INJECTION_UPDATE;
virtual ~WellState() = default;
@@ -69,11 +66,8 @@ public:
this->phase_usage_ = pu;
}
const WellMapType& wellMap() const { return wellMap_; }
WellMapType& wellMap() { return wellMap_; }
std::size_t size() const {
return this->wellMap_.size();
return this->wells_.size();
}
@@ -82,8 +76,6 @@ public:
return this->size();
}
int wellIndex(const std::string& wellName) const;
const ParallelWellInfo& parallelWellInfo(std::size_t well_index) const;
@@ -149,18 +141,6 @@ public:
/// One rate pr well
double brineWellRate(const int w) const;
const WellContainer<std::vector<double>>& wellReservoirRates() const { return well_reservoir_rates_; }
std::vector<double>& wellReservoirRates(std::size_t well_index)
{
return well_reservoir_rates_[well_index];
}
const std::vector<double>& wellReservoirRates(std::size_t well_index) const
{
return well_reservoir_rates_[well_index];
}
template<class Comm>
void communicateGroupRates(const Comm& comm);
@@ -255,32 +235,17 @@ public:
}
/// One rate per well and phase.
const WellContainer<std::vector<double>>& wellRates() const { return wellrates_; }
std::vector<double>& wellRates(std::size_t well_index) { return wellrates_[well_index]; }
const std::vector<double>& wellRates(std::size_t well_index) const { return wellrates_[well_index]; }
std::size_t numPerf(std::size_t well_index) const { return this->perfdata[well_index].size(); }
PerfData& perfData(const std::string& wname) {
return this->perfdata[wname];
}
const PerfData& perfData(const std::string& wname) const {
return this->perfdata[wname];
}
PerfData& perfData(std::size_t well_index) {
return this->perfdata[well_index];
}
const PerfData& perfData(std::size_t well_index) const {
return this->perfdata[well_index];
}
std::vector<double>& wellRates(std::size_t well_index) { return this->wells_[well_index].surface_rates; }
const std::vector<double>& wellRates(std::size_t well_index) const { return this->wells_[well_index].surface_rates; }
const std::string& name(std::size_t well_index) const {
return this->wells_.well_name(well_index);
}
std::optional<std::size_t> index(const std::string& well_name) const {
return this->wells_.well_index(well_name);
}
const SingleWellState& well(std::size_t well_index) const {
return this->wells_[well_index];
}
@@ -302,7 +267,6 @@ public:
}
private:
WellMapType wellMap_;
// Use of std::optional<> here is a technical crutch, the
// WellStateFullyImplicitBlackoil class should be default constructible,
// whereas the GlobalWellInfo is not.
@@ -313,18 +277,11 @@ private:
WellContainer<SingleWellState> wells_;
WellContainer<const ParallelWellInfo*> parallel_well_info_;
WellContainer<std::vector<double>> wellrates_;
WellContainer<PerfData> perfdata;
// The well_rates variable is defined for all wells on all processors. The
// bool in the value pair is whether the current process owns the well or
// not.
std::map<std::string, std::pair<bool, std::vector<double>>> well_rates;
// phase rates under reservoir condition for wells
// or voidage phase rates
WellContainer<std::vector<double>> well_reservoir_rates_;
data::Segment
reportSegmentResults(const PhaseUsage& pu,
const int well_id,
@@ -359,7 +316,6 @@ private:
const SummaryState& summary_state);
void initSingleWell(const std::vector<double>& cellPressures,
const int w,
const Well& well,
const std::vector<PerforationData>& well_perf_data,
const ParallelWellInfo* well_info,
tests/test_wellstate.cpp
+11 -19
View File
@@ -294,8 +294,8 @@ BOOST_AUTO_TEST_CASE(Pressure)
}
}
const auto& perf_data = wstate.perfData("PROD01");
const auto& ws = wstate.well("PROD01");
const auto& perf_data = ws.perf_data;
(void) perf_data;
}
@@ -362,7 +362,8 @@ BOOST_AUTO_TEST_CASE(STOP_well)
std::vector<Opm::ParallelWellInfo> pinfos;
auto wstate = buildWellState(setup, 0, pinfos);
for (std::size_t well_index = 0; well_index < setup.sched.numWells(0); well_index++) {
const auto& perf_data = wstate.perfData(well_index);
const auto& ws = wstate.well(well_index);
const auto& perf_data = ws.perf_data;
for (const auto& p : perf_data.pressure)
BOOST_CHECK(p > 0);
}
@@ -541,19 +542,10 @@ BOOST_AUTO_TEST_CASE(TESTSegmentState2) {
BOOST_AUTO_TEST_CASE(TESTPerfData) {
const auto& deck_string = R"(
RUNSPEC
OIL
WATER
GAS
)";
Opm::PhaseUsage pu = Opm::phaseUsageFromDeck(Opm::Parser{}.parseString(deck_string));
Opm::PerfData pd1(3, true, pu);
Opm::PerfData pd2(3, true, pu);
Opm::PerfData pd3(2, true, pu);
Opm::PerfData pd4(3, false, pu);
Opm::PerfData pd1(3, true, 3);
Opm::PerfData pd2(3, true, 3);
Opm::PerfData pd3(2, true, 3);
Opm::PerfData pd4(3, false, 3);
for (std::size_t i = 0; i < 3; i++) {
@@ -577,9 +569,9 @@ GAS
BOOST_AUTO_TEST_CASE(TestSingleWellState) {
Opm::SingleWellState ws1(true, 3, 1);
Opm::SingleWellState ws2(true, 3, 2);
Opm::SingleWellState ws3(false, 3, 3);
Opm::SingleWellState ws1(true, 10, 3, 1);
Opm::SingleWellState ws2(true, 10, 3, 2);
Opm::SingleWellState ws3(false, 10, 3, 3);
ws1.bhp = 100;
ws1.thp = 200;