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

Use well and group controls from opm-common.

Browse Source 
This PR remove the usage of well_control_ from opm-core
and instead uses the control classes for wells and groups
from opm-common.
This PR also removes the usage of the group classes from
opm-core.
This commit is contained in:
Tor Harald Sandve 2019-08-07 14:13:11 +02:00
parent ba4b965785
commit 53896ffca8
16 changed files with 2588 additions and 1569 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists_files.cmake
View File

@ -195,6 +195,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/VFPHelpers.hpp
opm/simulators/wells/VFPInjProperties.hpp
opm/simulators/wells/VFPProdProperties.hpp
opm/simulators/wells/WellGroupHelpers.hpp
opm/simulators/wells/WellHelpers.hpp
opm/simulators/wells/WellInterface.hpp
opm/simulators/wells/WellInterface_impl.hpp

6
opm/simulators/flow/BlackoilModelEbos.hpp
View File

@ -270,12 +270,6 @@ namespace Opm {
OPM_THROW(Opm::NumericalIssue, "Too large residual found!");
}
}
// checking whether the group targets are converged
if (wellModel().wellCollection().groupControlActive()) {
report.converged = report.converged && wellModel().wellCollection().groupTargetConverged(wellModel().wellState().wellRates());
}
report.update_time += perfTimer.stop();
residual_norms_history_.push_back(residual_norms);
if (!report.converged) {

1
opm/simulators/flow/MissingFeatures.cpp
View File

@ -422,7 +422,6 @@ namespace MissingFeatures {
"NOWARN",
"NSTACK",
"NUMRES",
"NUPCOL",
"OLDTRAN",
"OPERATER",
"OPTIONS",

32
opm/simulators/wells/BlackoilWellModel.hpp
View File

@ -47,6 +47,7 @@
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/simulators/wells/StandardWell.hpp>
#include <opm/simulators/wells/MultisegmentWell.hpp>
#include <opm/simulators/wells/WellGroupHelpers.hpp>
#include <opm/simulators/timestepping/gatherConvergenceReport.hpp>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
@ -203,11 +204,6 @@ namespace Opm {
// Check if well equations is converged.
ConvergenceReport getWellConvergence(const std::vector<Scalar>& B_avg) const;
// return all the wells.
const WellCollection& wellCollection() const;
// return non const reference to all the wells.
WellCollection& wellCollection();
// return the internal well state, ignore the passed one.
// Used by the legacy code to make it compatible with the legacy well models.
const WellState& wellState(const WellState& well_state OPM_UNUSED) const;
@ -306,9 +302,7 @@ namespace Opm {
// xw to update Well State
void recoverWellSolutionAndUpdateWellState(const BVector& x);
void updateWellControls(Opm::DeferredLogger& deferred_logger);
void updateGroupControls(Opm::DeferredLogger& deferred_logger);
void updateWellControls(Opm::DeferredLogger& deferred_logger, const bool checkGroupControl);
// setting the well_solutions_ based on well_state.
void updatePrimaryVariables(Opm::DeferredLogger& deferred_logger);
@ -320,17 +314,12 @@ namespace Opm {
void computeAverageFormationFactor(std::vector<Scalar>& B_avg) const;
void applyVREPGroupControl();
void computeWellVoidageRates(std::vector<double>& well_voidage_rates,
std::vector<double>& voidage_conversion_coeffs) const;
// Calculating well potentials for each well
void computeWellPotentials(std::vector<double>& well_potentials, Opm::DeferredLogger& deferred_logger);
void computeWellPotentials(std::vector<double>& well_potentials, const int reportStepIdx, Opm::DeferredLogger& deferred_logger);
const std::vector<double>& wellPerfEfficiencyFactors() const;
void calculateEfficiencyFactors();
void calculateEfficiencyFactors(const int reportStepIdx);
// it should be able to go to prepareTimeStep(), however, the updateWellControls() and initPrimaryVariablesEvaluation()
// makes it a little more difficult. unless we introduce if (iterationIdx != 0) to avoid doing the above functions
@ -350,18 +339,12 @@ namespace Opm {
int numPhases() const;
void resetWellControlFromState() const;
void assembleWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger);
// some preparation work, mostly related to group control and RESV,
// at the beginning of each time step (Not report step)
void prepareTimeStep(Opm::DeferredLogger& deferred_logger);
void prepareGroupControl(Opm::DeferredLogger& deferred_logger);
void computeRESV(Opm::DeferredLogger& deferred_logger);
void extractLegacyCellPvtRegionIndex_();
void extractLegacyDepth_();
@ -391,6 +374,13 @@ namespace Opm {
bool anyMSWellOpenLocal(const Wells* wells) const;
const Well2& getWellEcl(const std::string& well_name) const;
void checkGroupConstraints(const Group2& group, Opm::DeferredLogger& deferred_logger);
void actionOnBrokenConstraints(const Group2& group, const Group2::ExceedAction& exceed_action, const Group2::ProductionCMode& newControl, const int reportStepIdx, Opm::DeferredLogger& deferred_logger);
void actionOnBrokenConstraints(const Group2& group, const Group2::InjectionCMode& newControl, const int reportStepIdx, Opm::DeferredLogger& deferred_logger);
};

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

File diff suppressed because it is too large Load Diff

20
opm/simulators/wells/MultisegmentWell.hpp
View File

@ -124,7 +124,7 @@ namespace Opm
Opm::DeferredLogger& deferred_logger) const override;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const override;
virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const override;
/// Ax = Ax - C D^-1 B x
virtual void apply(const BVector& x, BVector& Ax) const override;
@ -188,7 +188,6 @@ namespace Opm
using Base::saturation_table_number_;
using Base::well_efficiency_factor_;
using Base::gravity_;
using Base::well_controls_;
using Base::perf_depth_;
using Base::num_components_;
using Base::connectionRates_;
@ -200,6 +199,8 @@ namespace Opm
using Base::ebosCompIdxToFlowCompIdx;
using Base::getAllowCrossFlow;
using Base::scalingFactor;
using Base::wellIsStopped_;
// TODO: trying to use the information from the Well opm-parser as much
// as possible, it will possibly be re-implemented later for efficiency reason.
@ -325,7 +326,7 @@ namespace Opm
EvalWell getSegmentRate(const int seg, const int comp_idx) const;
EvalWell getSegmentRateUpwinding(const int seg, const int comp_idx) const;
EvalWell getSegmentRateUpwinding(const int seg, const size_t comp_idx) const;
EvalWell getSegmentGTotal(const int seg) const;
@ -340,7 +341,10 @@ namespace Opm
std::vector<double>& well_flux,
Opm::DeferredLogger& deferred_logger);
void assembleControlEq(Opm::DeferredLogger& deferred_logger) const;
void assembleControlEq(const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, Opm::DeferredLogger& deferred_logger);
void assembleGroupProductionControl(const Group2& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
void assembleGroupInjectionControl(const Group2& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, const Well2::InjectorType& injectorType, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
void assemblePressureEq(const int seg) const;
@ -392,12 +396,14 @@ namespace Opm
void detectOscillations(const std::vector<double>& measure_history,
const int it, bool& oscillate, bool& stagnate) const;
double getResidualMeasureValue(const std::vector<double>& residuals,
double getResidualMeasureValue(const WellState& well_state,
const std::vector<double>& residuals,
DeferredLogger& deferred_logger) const;
double getControlTolerance(DeferredLogger& deferred_logger) const;
double getControlTolerance(const WellState& well_state, DeferredLogger& deferred_logger) const;
void checkConvergenceControlEq(ConvergenceReport& report,
void checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const;
void updateUpwindingSegments();

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

File diff suppressed because it is too large Load Diff

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

@ -163,7 +163,8 @@ namespace Opm
Opm::DeferredLogger& deferred_logger) const override;
/// check whether the well equations get converged for this well
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg,
virtual ConvergenceReport getWellConvergence(const WellState& well_state,
const std::vector<double>& B_avg,
Opm::DeferredLogger& deferred_logger) const override;
/// Ax = Ax - C D^-1 B x
@ -210,10 +211,9 @@ namespace Opm
using Base::wsolvent;
using Base::wpolymer;
using Base::wfoam;
using Base::wellHasTHPConstraints;
using Base::mostStrictBhpFromBhpLimits;
using Base::scalingFactor;
using Base::scaleProductivityIndex;
using Base::mostStrictBhpFromBhpLimits;
// protected member variables from the Base class
using Base::current_step_;
@ -232,7 +232,6 @@ namespace Opm
using Base::comp_frac_;
using Base::well_index_;
using Base::index_of_well_;
using Base::well_controls_;
using Base::well_type_;
using Base::num_components_;
using Base::connectionRates_;
@ -241,6 +240,8 @@ namespace Opm
using Base::perf_length_;
using Base::bore_diameters_;
using Base::wellIsStopped_;
// total number of the well equations and primary variables
// there might be extra equations be used, numWellEq will be updated during the initialization
int numWellEq_ = numStaticWellEq;
@ -364,7 +365,8 @@ namespace Opm
Opm::DeferredLogger& deferred_logger);
template <class ValueType>
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const int control_index, Opm::DeferredLogger& deferred_logger) const;
ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const Well2& well, const SummaryState& summaryState, Opm::DeferredLogger& deferred_logger) const;
double calculateThpFromBhp(const std::vector<double>& rates, const double bhp, Opm::DeferredLogger& deferred_logger) const;
@ -390,7 +392,11 @@ namespace Opm
void updateThp(WellState& well_state, Opm::DeferredLogger& deferred_logger) const;
void assembleControlEq(Opm::DeferredLogger& deferred_logger);
void assembleControlEq(const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, Opm::DeferredLogger& deferred_logger);
void assembleGroupProductionControl(const Group2& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
void assembleGroupInjectionControl(const Group2& group, const WellState& well_state, const Opm::Schedule& schedule, const SummaryState& summaryState, const Well2::InjectorType& injectorType, EvalWell& control_eq, double efficincyFactor, Opm::DeferredLogger& deferred_logger);
// handle the non reasonable fractions due to numerical overshoot
void processFractions() const;
@ -444,7 +450,7 @@ namespace Opm
// calculate the BHP from THP target based on IPR
// TODO: we need to check the operablility here first, if not operable, then maybe there is
// no point to do this
double calculateBHPWithTHPTargetIPR(Opm::DeferredLogger& deferred_logger) const;
double calculateBHPWithTHPTargetIPR(const Well2& well, const SummaryState& summaryState, Opm::DeferredLogger& deferred_logger) const;
// relaxation factor considering only one fraction value
static double relaxationFactorFraction(const double old_value,
@ -493,7 +499,8 @@ namespace Opm
virtual void updateWaterThroughput(const double dt, WellState& well_state) const override;
// checking the convergence of the well control equations
void checkConvergenceControlEq(ConvergenceReport& report,
void checkConvergenceControlEq(const WellState& well_state,
ConvergenceReport& report,
DeferredLogger& deferred_logger) const;
// checking convergence of extra equations, if there are any
@ -507,6 +514,7 @@ namespace Opm
const int perf,
DeferredLogger& deferred_logger);
};
}

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

File diff suppressed because it is too large Load Diff

216
opm/simulators/wells/WellGroupHelpers.hpp Normal file
View File

@ -0,0 +1,216 @@
/*
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_WELLGROUPHELPERS_HEADER_INCLUDED
#define OPM_WELLGROUPHELPERS_HEADER_INCLUDED
#include <vector>
namespace Opm {
namespace wellGroupHelpers
{
inline void setCmodeGroup(const Group2& group, const Schedule& schedule, const SummaryState& summaryState, const int reportStepIdx, WellStateFullyImplicitBlackoil& wellState) {
for (const std::string& groupName : group.groups()) {
setCmodeGroup( schedule.getGroup2(groupName, reportStepIdx), schedule, summaryState, reportStepIdx, wellState);
}
// use FLD as default control
wellState.setCurrentInjectionGroupControl(group.name(), Group2::InjectionCMode::FLD);
wellState.setCurrentProductionGroupControl(group.name(), Group2::ProductionCMode::FLD);
if (group.isInjectionGroup()) {
const auto controls = group.injectionControls(summaryState);
wellState.setCurrentInjectionGroupControl(group.name(), controls.cmode);
}
if (group.isProductionGroup()) {
const auto controls = group.productionControls(summaryState);
wellState.setCurrentProductionGroupControl(group.name(), controls.cmode);
}
}
inline void accumulateGroupEfficiencyFactor(const Group2& group, const Schedule& schedule, const int reportStepIdx, double& factor) {
factor *= group.getGroupEfficiencyFactor();
if (group.parent() != "FIELD")
accumulateGroupEfficiencyFactor(schedule.getGroup2(group.parent(), reportStepIdx), schedule, reportStepIdx, factor);
}
inline double sumWellRates(const Group2& group, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const int phasePos, const bool injector) {
double rate = 0.0;
for (const std::string& groupName : group.groups()) {
const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
rate += groupTmp.getGroupEfficiencyFactor()*sumWellRates(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];
const auto& wellEcl = schedule.getWell2(wellName, reportStepIdx);
//only count producers or injectors
if ( (wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector))
continue;
double factor = wellEcl.getEfficiencyFactor();
const auto wellrate_index = well_index * wellState.numPhases();
if (injector)
rate += factor * wellState.wellRates()[ wellrate_index + phasePos];
else
rate -= factor * wellState.wellRates()[ wellrate_index + phasePos];
}
return rate;
}
inline double sumWellResRates(const Group2& group, const Schedule& schedule, const WellStateFullyImplicitBlackoil& wellState, const int reportStepIdx, const int phasePos, const bool injector) {
double rate = 0.0;
for (const std::string& groupName : group.groups()) {
const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
rate += groupTmp.getGroupEfficiencyFactor() * sumWellResRates(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];
const auto& wellEcl = schedule.getWell2(wellName, reportStepIdx);
//only count producers or injectors
if ( (wellEcl.isProducer() && injector) || (wellEcl.isInjector() && !injector))
continue;
double factor = wellEcl.getEfficiencyFactor();
const auto wellrate_index = well_index * wellState.numPhases();
if (injector)
rate += factor * wellState.wellReservoirRates()[ wellrate_index + phasePos];
else
rate -= factor * wellState.wellReservoirRates()[ wellrate_index + phasePos];
}
return rate;
}
inline void setGroupControl(const Group2& group, const Schedule& schedule, const int reportStepIdx, const bool injector, WellStateFullyImplicitBlackoil& wellState) {
for (const std::string& groupName : group.groups()) {
const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
setGroupControl(groupTmp, schedule, reportStepIdx, injector, wellState);
if (injector)
wellState.setCurrentInjectionGroupControl(groupName, Group2::InjectionCMode::FLD);
else
wellState.setCurrentProductionGroupControl(groupName, Group2::ProductionCMode::FLD);
}
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];
const auto& wellEcl = schedule.getWell2(wellName, reportStepIdx);
if (!wellEcl.isAvailableForGroupControl())
continue;
if (wellEcl.isProducer() && !injector)
wellState.currentProductionControls()[well_index] = Well2::ProducerCMode::GRUP;
if (wellEcl.isInjector() && injector)
wellState.currentInjectionControls()[well_index] = Well2::InjectorCMode::GRUP;
}
}
inline double computeGuideRate(const WellStateFullyImplicitBlackoil& wellState, const Well2& well, const int phasePos )
{
const auto& end = wellState.wellMap().end();
const auto& it = wellState.wellMap().find( well.name() );
if (it == end) // the well is not found
return 0.0;
int well_index = it->second[0];
Opm::Well2::GuideRateTarget GuideRatePhase = well.getGuideRatePhase();
#warning handle the case when phase != GuideRatePhase
if (GuideRatePhase == Opm::Well2::GuideRateTarget::UNDEFINED) {
const auto wellrate_index = well_index * wellState.numPhases();
double guideRate = wellState.wellPotentials()[wellrate_index + phasePos];
//if (guideRate == 0)
// guideRate = 1; // if the well potential is not defined set it to 1
#warning Add support for GUIDERAT
return guideRate*well.getGuideRateScalingFactor();
}
return well.getGuideRate()*well.getGuideRateScalingFactor();
}
inline void computeTotalGuideRate(const Group2& group, const WellStateFullyImplicitBlackoil& wellState, const Schedule& schedule, const int reportStepIdx, const int phasePos, const bool isInjector, double& groupTotalGuideRate, double& groupTargetReduction )
{
for (const std::string& groupName : group.groups()) {
const Group2& groupTmp = schedule.getGroup2(groupName, reportStepIdx);
computeTotalGuideRate(groupTmp, wellState, schedule, reportStepIdx, phasePos, isInjector, groupTotalGuideRate, groupTargetReduction);
}
#warning return groupGuideRate if set, else sum it from below
for (const std::string& wellName : group.wells()) {
const auto& wellTmp = schedule.getWell2(wellName, reportStepIdx);
if (wellTmp.isProducer() && isInjector)
continue;
if (wellTmp.isInjector() && !isInjector)
continue;
const auto& end = wellState.wellMap().end();
const auto& it = wellState.wellMap().find( wellName );
if (it == end) // the well is not found
continue;
int well_index = it->second[0];
const auto wellrate_index = well_index * wellState.numPhases();
// only include wells under group control
if (isInjector) {
if (wellState.currentInjectionControls()[well_index] == Well2::InjectorCMode::GRUP)
groupTotalGuideRate += computeGuideRate(wellState, wellTmp, phasePos);
else
groupTargetReduction += wellState.wellRates()[wellrate_index + phasePos];
} else {
if (wellState.currentProductionControls()[well_index] == Well2::ProducerCMode::GRUP)
groupTotalGuideRate += computeGuideRate(wellState, wellTmp, phasePos);
else
groupTargetReduction -= wellState.wellRates()[wellrate_index + phasePos];
}
}
}
} // namespace wellGroupHelpers
}
#endif

85
opm/simulators/wells/WellHelpers.hpp
View File

@ -36,90 +36,7 @@ namespace Opm {
{
inline
double rateToCompare(const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
const double* distr)
{
double rate = 0.0;
for (int phase = 0; phase < num_phases; ++phase) {
// Important: well_phase_flow_rate is ordered with all phase rates for first
// well first, then all phase rates for second well etc.
rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase];
}
return rate;
}
inline
bool constraintBroken(const std::vector<double>& bhp,
const std::vector<double>& thp,
const std::vector<double>& well_phase_flow_rate,
const int well,
const int num_phases,
const WellType& well_type,
const WellControls* wc,
const int ctrl_index)
{
const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index);
const double target = well_controls_iget_target(wc, ctrl_index);
const double* distr = well_controls_iget_distr(wc, ctrl_index);
bool broken = false;
switch (well_type) {
case INJECTOR:
{
switch (ctrl_type) {
case BHP:
broken = bhp[well] > target;
break;
case THP:
broken = thp[well] > target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) > target;
break;
}
}
break;
case PRODUCER:
{
switch (ctrl_type) {
case BHP:
broken = bhp[well] < target;
break;
case THP:
broken = thp[well] < target;
break;
case RESERVOIR_RATE: // Intentional fall-through
case SURFACE_RATE:
// Note that the rates compared below are negative,
// so breaking the constraints means: too high flow rate
// (as for injection).
broken = rateToCompare(well_phase_flow_rate,
well, num_phases, distr) < target;
break;
}
}
break;
default:
OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells.");
}
return broken;
}
// --------- Types ---------
// --------- Types ---------
/**
* Simple hydrostatic correction for VFP table

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

@ -2,6 +2,7 @@
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).
@ -39,6 +40,7 @@
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/VFPProperties.hpp>
#include <opm/simulators/wells/WellHelpers.hpp>
#include <opm/simulators/wells/WellGroupHelpers.hpp>
#include <opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp>
#include <opm/simulators/flow/BlackoilModelParametersEbos.hpp>
@ -148,7 +150,7 @@ namespace Opm
virtual void initPrimaryVariablesEvaluation() const = 0;
virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const = 0;
virtual ConvergenceReport getWellConvergence(const WellState& well_state, const std::vector<double>& B_avg, Opm::DeferredLogger& deferred_logger) const = 0;
virtual void solveEqAndUpdateWellState(WellState& well_state, Opm::DeferredLogger& deferred_logger) = 0;
@ -232,7 +234,7 @@ namespace Opm
void closeCompletions(WellTestState& wellTestState);
const Well2* wellEcl() const;
const Well2& wellEcl() const;
// TODO: theoretically, it should be a const function
// Simulator is not const is because that assembleWellEq is non-const Simulator
@ -250,7 +252,7 @@ namespace Opm
bool isOperable() const;
/// Returns true if the well has one or more THP limits/constraints.
bool wellHasTHPConstraints() const;
bool wellHasTHPConstraints(const SummaryState& summaryState) const;
/// Returns true if the well is currently in prediction mode (i.e. not history mode).
bool underPredictionMode() const;
@ -258,6 +260,17 @@ namespace Opm
// update perforation water throughput based on solved water rate
virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0;
void stopWell() {
wellIsStopped_ = true;
}
void openWell() {
wellIsStopped_ = false;
}
bool wellIsStopped() {
return wellIsStopped_;
}
protected:
// to indicate a invalid completion
@ -284,9 +297,6 @@ namespace Opm
// typically, it should apply to injection wells
std::vector<double> comp_frac_;
// controls for this well
struct WellControls* well_controls_;
// number of the perforations for this well
int number_of_perforations_;
@ -356,6 +366,8 @@ namespace Opm
std::vector<RateVector> connectionRates_;
bool wellIsStopped_;
const PhaseUsage& phaseUsage() const;
int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
@ -372,14 +384,14 @@ namespace Opm
const WellState& well_state,
Opm::DeferredLogger& deferred_logger) const;
double getTHPConstraint(Opm::DeferredLogger& deferred_logger) const;
double getTHPConstraint(const SummaryState& summaryState) const;
int getControlIndex(const WellControlType& type) const;
// Component fractions for each phase for the well
const std::vector<double>& compFrac() const;
double mostStrictBhpFromBhpLimits(Opm::DeferredLogger& deferred_logger) const;
double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
struct RatioLimitCheckReport;
@ -437,9 +449,9 @@ namespace Opm
WellTestState& well_test_state,
Opm::DeferredLogger& deferred_logger) const;
void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state,
const std::vector<double>& B_avg,
Opm::DeferredLogger& deferred_logger);
void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state,
const std::vector<double>& B_avg,
Opm::DeferredLogger& deferred_logger);
bool solveWellEqUntilConverged(const Simulator& ebosSimulator,
const std::vector<double>& B_avg,
@ -450,18 +462,18 @@ namespace Opm
void initCompletions();
WellControls* createWellControlsWithBHPAndTHP(DeferredLogger& deferred_logger) const;
// count the number of times an output log message is created in the productivity
// index calculations
int well_productivity_index_logger_counter_;
bool checkConstraints(WellState& well_state, const SummaryState& summaryState);
};
// definition of the struct OperabilityStatus
template<typename TypeTag>
struct

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

@ -73,8 +73,6 @@ namespace Opm
wells->comp_frac + index_begin + number_of_phases_, comp_frac_.begin() );
}
well_controls_ = wells->ctrls[index_well];
ref_depth_ = wells->depth_ref[index_well];
// perforations related
@ -109,6 +107,12 @@ namespace Opm
well_productivity_index_logger_counter_ = 0;
wellIsStopped_ = false;
if (well.getStatus() == Well2::Status::STOP) {
wellIsStopped_ = true;
}
}
template<typename TypeTag>
@ -199,14 +203,6 @@ namespace Opm
template<typename TypeTag>
WellControls*
WellInterface<TypeTag>::
wellControls() const
{
return well_controls_;
}
template<typename TypeTag>
int
WellInterface<TypeTag>::
@ -242,11 +238,11 @@ namespace Opm
template<typename TypeTag>
const Well2*
const Well2&
WellInterface<TypeTag>::
wellEcl() const
{
return std::addressof(well_ecl_);
return well_ecl_;
}
@ -372,107 +368,69 @@ namespace Opm
template<typename TypeTag>
double
WellInterface<TypeTag>::
mostStrictBhpFromBhpLimits(Opm::DeferredLogger& deferred_logger) const
{
double bhp;
// initial bhp value, making the value not usable
switch( well_type_ ) {
case INJECTOR:
bhp = std::numeric_limits<double>::max();
break;
case PRODUCER:
bhp = -std::numeric_limits<double>::max();
break;
default:
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name(), deferred_logger);
}
// The number of the well controls/constraints
const int nwc = well_controls_get_num(well_controls_);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
// finding a BHP constraint
if (well_controls_iget_type(well_controls_, ctrl_index) == BHP) {
// get the bhp constraint value, it should always be postive assummingly
const double bhp_target = well_controls_iget_target(well_controls_, ctrl_index);
switch(well_type_) {
case INJECTOR: // using the lower bhp contraint from Injectors
if (bhp_target < bhp) {
bhp = bhp_target;
}
break;
case PRODUCER:
if (bhp_target > bhp) {
bhp = bhp_target;
}
break;
default:
OPM_DEFLOG_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type for well " << name(), deferred_logger);
} // end of switch
}
}
return bhp;
}
template<typename TypeTag>
bool
WellInterface<TypeTag>::
wellHasTHPConstraints() const
wellHasTHPConstraints(const SummaryState& summaryState) const
{
return getControlIndex(THP) >= 0;
if (well_ecl_.isInjector()) {
const auto controls = well_ecl_.injectionControls(summaryState);
if (controls.hasControl(Well2::InjectorCMode::THP))
return true;
}
if (well_ecl_.isProducer( )) {
const auto controls = well_ecl_.productionControls(summaryState);
if (controls.hasControl(Well2::ProducerCMode::THP))
return true;
}
return false;
}
template<typename TypeTag>
double
WellInterface<TypeTag>::
getTHPConstraint(Opm::DeferredLogger& deferred_logger) const
mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
{
const int thp_control_index = getControlIndex(THP);
if (thp_control_index < 0) {
OPM_DEFLOG_THROW(std::runtime_error, " there is no THP constraint/limit for well " << name()
<< ", while we are requesting it ", deferred_logger);
if (well_ecl_.isInjector()) {
const auto& controls = well_ecl_.injectionControls(summaryState);
return controls.bhp_limit;
}
return well_controls_iget_target(well_controls_, thp_control_index);
if (well_ecl_.isProducer( )) {
const auto& controls = well_ecl_.productionControls(summaryState);
return controls.bhp_limit;
}
return 0.0;
}
template<typename TypeTag>
int
double
WellInterface<TypeTag>::
getControlIndex(const WellControlType& type) const
getTHPConstraint(const SummaryState& summaryState) const
{
const int nwc = well_controls_get_num(well_controls_);
for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(well_controls_, ctrl_index) == type) {
return ctrl_index;
}
if (well_ecl_.isInjector()) {
const auto& controls = well_ecl_.injectionControls(summaryState);
return controls.thp_limit;
}
return -1;
if (well_ecl_.isProducer( )) {
const auto& controls = well_ecl_.productionControls(summaryState);
return controls.thp_limit;
}
return 0.0;
}
template<typename TypeTag>
void
WellInterface<TypeTag>::
@ -480,81 +438,36 @@ namespace Opm
WellState& well_state,
Opm::DeferredLogger& deferred_logger) /* const */
{
const auto& summaryState = ebos_simulator.vanguard().summaryState();
const auto& well = well_ecl_;
std::string from;
if (well.isInjector()) {
from = Well2::InjectorCMode2String(well_state.currentInjectionControls()[index_of_well_]);
} else {
from = Well2::ProducerCMode2String(well_state.currentProductionControls()[index_of_well_]);
}
bool changed = checkConstraints(well_state, summaryState);
auto cc = Dune::MPIHelper::getCollectiveCommunication();
const int np = number_of_phases_;
const int w = index_of_well_;
const int old_control_index = well_state.currentControls()[w];
// Find, for each well, if any constraints are broken. If so,
// switch control to first broken constraint.
WellControls* wc = well_controls_;
// Loop over all controls except the current one, and also
// skip any RESERVOIR_RATE controls, since we cannot
// handle those.
const int nwc = well_controls_get_num(wc);
// the current control index
int current = well_state.currentControls()[w];
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (ctrl_index == current) {
// This is the currently used control, so it is
// used as an equation. So this is not used as an
// inequality constraint, and therefore skipped.
continue;
}
if (wellhelpers::constraintBroken(
well_state.bhp(), well_state.thp(), well_state.wellRates(),
w, np, well_type_, wc, ctrl_index)) {
// if the well can not work under THP / BHP control, we should not switch to THP / BHP control
const bool cannot_switch_to_bhp = well_controls_iget_type(wc, ctrl_index) == BHP && !operability_status_.isOperableUnderBHPLimit();
const bool cannot_switch_to_thp = well_controls_iget_type(wc, ctrl_index) == THP && !operability_status_.isOperableUnderTHPLimit();
const bool cannot_switch = cannot_switch_to_bhp || cannot_switch_to_thp;
if ( !cannot_switch ) {
// ctrl_index will be the index of the broken constraint after the loop.
break;
} else {
// before we figure out to handle it, we give some debug information here
if ( well_controls_iget_type(wc, ctrl_index) == BHP && !operability_status_.isOperableUnderBHPLimit() ) {
deferred_logger.debug("well " + name() + " breaks the BHP limit, while it is not operable under BHP limit");
}
if ( well_controls_iget_type(wc, ctrl_index) == THP && !operability_status_.isOperableUnderTHPLimit() ) {
deferred_logger.debug("well " + name() + " breaks the THP limit, while it is not operable under THP limit");
}
}
}
}
if (ctrl_index != nwc) {
// Constraint number ctrl_index was broken, switch to it.
well_state.currentControls()[w] = ctrl_index;
current = well_state.currentControls()[w];
well_controls_set_current( wc, current);
}
// the new well control indices after all the related updates,
const int updated_control_index = well_state.currentControls()[w];
// checking whether control changed
if (updated_control_index != old_control_index) {
auto from = well_controls_iget_type(wc, old_control_index);
auto to = well_controls_iget_type(wc, updated_control_index);
if (changed) {
std::string to;
if (well.isInjector()) {
to = Well2::InjectorCMode2String(well_state.currentInjectionControls()[index_of_well_]);
} else {
to = Well2::ProducerCMode2String(well_state.currentProductionControls()[index_of_well_]);
}
std::ostringstream ss;
ss << " Switching control mode for well " << name()
<< " from " << modestring[from]
<< " to " << modestring[to];
<< " from " << from
<< " to " << to;
if (cc.size() > 1) {
ss << " on rank " << cc.rank();
}
deferred_logger.info(ss.str());
}
if (updated_control_index != old_control_index) { // || well_collection_->groupControlActive()) {
updateWellStateWithTarget(ebos_simulator, well_state, deferred_logger);
updatePrimaryVariables(well_state, deferred_logger);
}
@ -886,7 +799,7 @@ namespace Opm
WellTestState& well_test_state,
Opm::DeferredLogger& deferred_logger) const
{
if (!isOperable()) {
if (!isOperable() && !wellIsStopped_) {
well_test_state.closeWell(name(), WellTestConfig::Reason::PHYSICAL, simulation_time);
if (write_message_to_opmlog) {
// TODO: considering auto shut in?
@ -911,6 +824,9 @@ namespace Opm
WellTestState& well_test_state,
Opm::DeferredLogger& deferred_logger) const
{
if (wellIsStopped_)
return;
const WellEconProductionLimits& econ_production_limits = well_ecl_.getEconLimits();
// if no limit is effective here, then continue to the next well
@ -1005,7 +921,7 @@ namespace Opm
if (write_message_to_opmlog) {
if (well_ecl_.getAutomaticShutIn()) {
const std::string msg = name() + std::string(" will be shut due to last completion closed");
deferred_logger.info(msg);
deferred_logger.info(msg);
} else {
const std::string msg = name() + std::string(" will be stopped due to last completion closed");
deferred_logger.info(msg);
@ -1204,19 +1120,6 @@ namespace Opm
double
WellInterface<TypeTag>::scalingFactor(const int phaseIdx) const
{
const WellControls* wc = well_controls_;
if (well_controls_get_current(wc) != -1 && well_controls_get_current_type(wc) == RESERVOIR_RATE) {
if (has_solvent && phaseIdx == contiSolventEqIdx ) {
typedef Opm::BlackOilSolventModule<TypeTag> SolventModule;
double coeff = 0;
rateConverter_.template calcCoeffSolvent<SolventModule>(0, pvtRegionIdx_, coeff);
return coeff;
}
// TODO: use the rateConverter here as well.
const double* distr = well_controls_get_current_distr(wc);
return distr[phaseIdx];
}
const auto& pu = phaseUsage();
if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) && pu.phase_pos[Water] == phaseIdx)
return 1.0;
@ -1292,7 +1195,7 @@ namespace Opm
do {
assembleWellEq(ebosSimulator, B_avg, dt, well_state, deferred_logger);
auto report = getWellConvergence(B_avg, deferred_logger);
auto report = getWellConvergence(well_state, B_avg, deferred_logger);
converged = report.converged();
if (converged) {
@ -1302,7 +1205,8 @@ namespace Opm
++it;
solveEqAndUpdateWellState(well_state, deferred_logger);
updateWellControl(ebosSimulator, well_state, deferred_logger);
// We don't allow for switching well controls while computing well potentials and testing wells
// updateWellControl(ebosSimulator, well_state, deferred_logger);
initPrimaryVariablesEvaluation();
} while (it < max_iter);
@ -1441,35 +1345,197 @@ namespace Opm
return operability_status_.isOperable();
}
template<typename TypeTag>
WellControls*
bool
WellInterface<TypeTag>::
createWellControlsWithBHPAndTHP(DeferredLogger& deferred_logger) const
{
WellControls* wc = well_controls_create();
well_controls_assert_number_of_phases(wc, number_of_phases_);
checkConstraints(WellState& well_state, const SummaryState& summaryState) {
// a well always has a bhp limit
const double invalid_alq = -1e100;
const double invalid_vfp = -2147483647;
const double bhp_limit = this->mostStrictBhpFromBhpLimits(deferred_logger);
well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc);
const auto& well = well_ecl_;
const PhaseUsage& pu = phaseUsage();
const int well_index = index_of_well_;
const auto wellrate_index = well_index * pu.num_phases;
bool changed = false;
if (this->wellHasTHPConstraints()) {
// it might be better to do it through EclipseState?
const double thp_limit = this->getTHPConstraint(deferred_logger);
const double thp_control_index = this->getControlIndex(THP);
const int vfp_number = well_controls_iget_vfp(well_controls_, thp_control_index);
const double alq = well_controls_iget_alq(well_controls_, thp_control_index);
well_controls_add_new(THP, thp_limit, alq, vfp_number, NULL, wc);
// // Stopped wells can not change control
// if (currentControl == "STOP")
// return newControl;
if (well.isInjector()) {
const auto controls = well.injectionControls(summaryState);
Opm::Well2::InjectorCMode& currentControl = well_state.currentInjectionControls()[well_index];
if (controls.hasControl(Well2::InjectorCMode::RATE) && currentControl != Well2::InjectorCMode::RATE)
{
Well2::InjectorType injectorType = controls.injector_type;
double current_rate = 0.0;
switch (injectorType) {
case Well2::InjectorType::WATER:
{
current_rate = well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Aqua] ];
break;
}
case Well2::InjectorType::OIL:
{
current_rate = well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Liquid] ];
break;
}
case Well2::InjectorType::GAS:
{
current_rate = well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ];
break;
}
default:
throw("Expected WATER, OIL or GAS as type for injectors " + well.name());
}
if (controls.surface_rate < current_rate) {
currentControl = Well2::InjectorCMode::RATE;
changed = true;
}
}
if (controls.hasControl(Well2::InjectorCMode::RESV) && currentControl != Well2::InjectorCMode::RESV)
{
double current_rate = 0.0;
if( pu.phase_used[BlackoilPhases::Aqua] )
current_rate += well_state.wellReservoirRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Aqua] ];
if( pu.phase_used[BlackoilPhases::Liquid] )
current_rate += well_state.wellReservoirRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Liquid] ];
if( pu.phase_used[BlackoilPhases::Vapour] )
current_rate += well_state.wellReservoirRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.reservoir_rate < current_rate) {
currentControl = Well2::InjectorCMode::RESV;
changed = true;
}
}
if (controls.hasControl(Well2::InjectorCMode::BHP) && currentControl != Well2::InjectorCMode::BHP)
{
const auto& bhp = controls.bhp_limit;
double current_bhp = well_state.bhp()[well_index];
if (bhp < current_bhp) {
currentControl = Well2::InjectorCMode::BHP;
changed = true;
}
}
if (controls.hasControl(Well2::InjectorCMode::THP) && currentControl != Well2::InjectorCMode::THP)
{
const auto& thp = controls.thp_limit;
double current_thp = well_state.thp()[well_index];
if (thp < current_thp) {
currentControl = Well2::InjectorCMode::THP;
changed = true;
}
}
}
return wc;
if (well.isProducer( )) {
const auto controls = well.productionControls(summaryState);
Well2::ProducerCMode& currentControl = well_state.currentProductionControls()[well_index];
if (controls.hasControl(Well2::ProducerCMode::ORAT) && currentControl != Well2::ProducerCMode::ORAT) {
double current_rate = -well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Liquid] ];
if (controls.oil_rate < current_rate ) {
currentControl = Well2::ProducerCMode::ORAT;
changed = true;
}
}
if (controls.hasControl(Well2::ProducerCMode::WRAT) && currentControl != Well2::ProducerCMode::WRAT ) {
double current_rate = -well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Aqua] ];
if (controls.water_rate < current_rate ) {
currentControl = Well2::ProducerCMode::WRAT;
changed = true;
}
}
if (controls.hasControl(Well2::ProducerCMode::GRAT) && currentControl != Well2::ProducerCMode::GRAT ) {
double current_rate = -well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.gas_rate < current_rate ) {
currentControl = Well2::ProducerCMode::GRAT;
changed = true;
}
}
if (controls.hasControl(Well2::ProducerCMode::LRAT) && currentControl != Well2::ProducerCMode::LRAT) {
double current_rate = -well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Liquid] ];
current_rate -= well_state.wellRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Aqua] ];
if (controls.liquid_rate < current_rate ) {
currentControl = Well2::ProducerCMode::LRAT;
changed = true;
}
}
if (controls.hasControl(Well2::ProducerCMode::RESV) && currentControl != Well2::ProducerCMode::RESV ) {
double current_rate = 0.0;
if( pu.phase_used[BlackoilPhases::Aqua] )
current_rate -= well_state.wellReservoirRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Aqua] ];
if( pu.phase_used[BlackoilPhases::Liquid] )
current_rate -= well_state.wellReservoirRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Liquid] ];
if( pu.phase_used[BlackoilPhases::Vapour] )
current_rate -= well_state.wellReservoirRates()[ wellrate_index + pu.phase_pos[BlackoilPhases::Vapour] ];
if (controls.prediction_mode && controls.resv_rate > current_rate) {
currentControl = Well2::ProducerCMode::RESV;
changed = true;
}
if (!controls.prediction_mode) {
const int fipreg = 0; // not considering the region for now
const int np = number_of_phases_;
std::vector<double> surface_rates(np, 0.0);
if( pu.phase_used[BlackoilPhases::Aqua] )
surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = controls.water_rate;
if( pu.phase_used[BlackoilPhases::Liquid] )
surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = controls.oil_rate;
if( pu.phase_used[BlackoilPhases::Vapour] )
surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = controls.gas_rate;
std::vector<double> voidage_rates(np, 0.0);
rateConverter_.calcReservoirVoidageRates(fipreg, pvtRegionIdx_, surface_rates, voidage_rates);
double resv_rate = 0.0;
for (int p = 0; p < np; ++p) {
resv_rate += voidage_rates[p];
}
if (resv_rate < current_rate) {
currentControl = Well2::ProducerCMode::RESV;
changed = true;
}
}
}
if (controls.hasControl(Well2::ProducerCMode::BHP) && currentControl != Well2::ProducerCMode::BHP )
{
const double bhp = controls.bhp_limit;
double current_bhp = well_state.bhp()[well_index];
if (bhp > current_bhp) {
currentControl = Well2::ProducerCMode::BHP;
changed = true;
}
}
if (controls.hasControl(Well2::ProducerCMode::THP) && currentControl != Well2::ProducerCMode::THP)
{
const auto& thp = controls.thp_limit;
double current_thp = well_state.thp()[well_index];
if (thp > current_thp) {
currentControl = Well2::ProducerCMode::THP;
changed = true;
}
}
}
return changed;
}

260
opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
View File

@ -38,6 +38,7 @@
#include <map>
#include <algorithm>
#include <array>
#include <iostream>
namespace Opm
{
@ -157,11 +158,12 @@ namespace Opm
}
}
current_controls_.resize(nw);
// The controls set in the Wells (specified in the DECK) are treated as default initial value
for (int w = 0; w < nw; ++w) {
current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
}
current_injection_controls_.resize(nw);
current_production_controls_.resize(nw);
current_injection_group_controls_.clear();
current_production_group_controls_.clear();
perfRateSolvent_.clear();
perfRateSolvent_.resize(nperf, 0.0);
productivity_index_.resize(nw * np, 0.0);
@ -186,12 +188,16 @@ namespace Opm
// thp
thp()[ newIndex ] = prevState->thp()[ oldIndex ];
// if there is no effective control event happens to the well, we use the current_controls_ from prevState
// Currently this is taken care of by updateWellStateFromTarge. Maybe we should just remove the initialization and just use updateWellStateFromTarget
//if (effective_events_occurred_[w]) {
// continue;
//}
// if there is no effective control event happens to the well, we use the current_injection/production_controls_ from prevState
// otherwise, we use the control specified in the deck
if (!effective_events_occurred_[w]) {
current_controls_[ newIndex ] = prevState->currentControls()[ oldIndex ];
// also change the one in the WellControls
well_controls_set_current(wells->ctrls[w], current_controls_[ newIndex ]);
current_injection_controls_[ newIndex ] = prevState->currentInjectionControls()[ oldIndex ];
current_production_controls_[ newIndex ] = prevState->currentProductionControls()[ oldIndex ];
}
// wellrates
@ -307,196 +313,44 @@ namespace Opm
}
}
/// Allocate and initialize if wells is non-null. Also tries
/// to give useful initial values to the bhp(), wellRates()
/// and perfPhaseRates() fields, depending on controls.
///
/// this method is only for flow_legacy!
template <class PrevWellState>
void initLegacy(const Wells* wells, const std::vector<double>& cellPressures , const PrevWellState& prevState, const PhaseUsage& pu)
{
// call init on base class
BaseType :: init(wells, cellPressures);
// if there are no well, do nothing in init
if (wells == 0) {
return;
}
const int nw = wells->number_of_wells;
if( nw == 0 ) return ;
// Initialize perfphaserates_, which must be done here.
const int np = wells->number_of_phases;
const int nperf = wells->well_connpos[nw];
well_reservoir_rates_.resize(nw * np, 0.0);
well_dissolved_gas_rates_.resize(nw, 0.0);
well_vaporized_oil_rates_.resize(nw, 0.0);
productivity_index_.resize(nw * np, 0.0);
well_potentials_.resize(nw*np, 0.0);
// Ensure that we start out with zero rates by default.
perfphaserates_.clear();
perfphaserates_.resize(nperf * np, 0.0);
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
const WellControls* ctrl = wells->ctrls[w];
if (well_controls_well_is_stopped(ctrl)) {
// Shut well: perfphaserates_ are all zero.
} else {
const int num_perf_this_well = wells->well_connpos[w + 1] - wells->well_connpos[w];
// Open well: Initialize perfphaserates_ to well
// rates divided by the number of perforations.
for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
for (int p = 0; p < np; ++p) {
perfphaserates_[np*perf + p] = wellRates()[np*w + p] / double(num_perf_this_well);
}
perfPress()[perf] = cellPressures[wells->well_cells[perf]];
}
}
}
// Initialize current_controls_.
// The controls set in the Wells object are treated as defaults,
// and also used for initial values.
current_controls_.resize(nw);
for (int w = 0; w < nw; ++w) {
current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
}
perfRateSolvent_.clear();
perfRateSolvent_.resize(nperf, 0.0);
// intialize wells that have been there before
// order may change so the mapping is based on the well name
if( ! prevState.wellMap().empty() )
{
typedef typename WellMapType :: const_iterator const_iterator;
const_iterator end = prevState.wellMap().end();
for (int w = 0; w < nw; ++w) {
std::string name( wells->name[ w ] );
const_iterator it = prevState.wellMap().find( name );
if( it != end )
{
const int oldIndex = (*it).second[ 0 ];
const int newIndex = w;
// bhp
bhp()[ newIndex ] = prevState.bhp()[ oldIndex ];
// thp
thp()[ newIndex ] = prevState.thp()[ oldIndex ];
// wellrates
for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
{
wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
}
// perfPhaseRates
const int oldPerf_idx_beg = (*it).second[ 1 ];
const int num_perf_old_well = (*it).second[ 2 ];
const int num_perf_this_well = wells->well_connpos[newIndex + 1] - wells->well_connpos[newIndex];
// copy perforation rates when the number of perforations is equal,
// otherwise initialize perfphaserates to well rates divided by the number of perforations.
if( num_perf_old_well == num_perf_this_well )
{
int old_perf_phase_idx = oldPerf_idx_beg *np;
for (int perf_phase_idx = wells->well_connpos[ newIndex ]*np;
perf_phase_idx < wells->well_connpos[ newIndex + 1]*np; ++perf_phase_idx, ++old_perf_phase_idx )
{
perfPhaseRates()[ perf_phase_idx ] = prevState.perfPhaseRates()[ old_perf_phase_idx ];
}
} else {
for (int perf = wells->well_connpos[newIndex]; perf < wells->well_connpos[newIndex + 1]; ++perf) {
for (int p = 0; p < np; ++p) {
perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(num_perf_this_well);
}
}
}
// perfPressures
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{
perfPress()[ perf ] = prevState.perfPress()[ oldPerf_idx ];
}
}
// perfSolventRates
if (pu.has_solvent) {
if( num_perf_old_well == num_perf_this_well )
{
int oldPerf_idx = oldPerf_idx_beg;
for (int perf = wells->well_connpos[ newIndex ];
perf < wells->well_connpos[ newIndex + 1]; ++perf, ++oldPerf_idx )
{
perfRateSolvent()[ perf ] = prevState.perfRateSolvent()[ oldPerf_idx ];
}
}
}
}
// If in the new step, there is no THP related target/limit anymore, its thp value should be
// set to zero.
const WellControls* ctrl = wells->ctrls[w];
const int nwc = well_controls_get_num(ctrl);
int ctrl_index = 0;
for (; ctrl_index < nwc; ++ctrl_index) {
if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
break;
}
}
// not finding any thp related control/limits
if (ctrl_index == nwc) {
thp()[w] = 0.;
}
}
}
{
// we need to create a trival segment related values to avoid there will be some
// multi-segment wells added later.
nseg_ = nw;
top_segment_index_.resize(nw);
seg_number_.resize(nw);
for (int w = 0; w < nw; ++w) {
top_segment_index_[w] = w;
seg_number_[w] = 1; // Top segment is segment #1
}
segpress_ = bhp();
segrates_ = wellRates();
}
}
// this method is only for flow_legacy!
template <class State, class PrevWellState>
void initLegacy(const Wells* wells, const State& state, const PrevWellState& prevState, const PhaseUsage& pu)
{
initLegacy(wells, state.pressure(), prevState, pu);
}
// this method is only for flow_legacy!
template <class State>
void resizeLegacy(const Wells* wells, const State& state, const PhaseUsage& pu)
{
const WellStateFullyImplicitBlackoil dummy_state{}; // Init with an empty previous state only resizes
initLegacy(wells, state, dummy_state, pu) ;
}
/// One rate per phase and well connection.
std::vector<double>& perfPhaseRates() { return perfphaserates_; }
const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
/// One current control per well.
std::vector<int>& currentControls() { return current_controls_; }
const std::vector<int>& currentControls() const { return current_controls_; }
/// One current control per injecting well.
std::vector<Opm::Well2::InjectorCMode>& currentInjectionControls() { return current_injection_controls_; }
const std::vector<Opm::Well2::InjectorCMode>& currentInjectionControls() const { return current_injection_controls_; }
/// One current control per producing well.
std::vector<Well2::ProducerCMode>& currentProductionControls() { return current_production_controls_; }
const std::vector<Well2::ProducerCMode>& currentProductionControls() const { return current_production_controls_; }
/// One current control per group.
void setCurrentProductionGroupControl(const std::string& groupName, const Group2::ProductionCMode& groupControl ) {
current_production_group_controls_[groupName] = groupControl;
}
const Group2::ProductionCMode& currentProductionGroupControl(const std::string& groupName) const {
auto it = current_production_group_controls_.find(groupName);
if (it == current_production_group_controls_.end())
OPM_THROW(std::logic_error, "Could not find any well control for " << groupName);
return it->second;
}
/// One current control per group.
void setCurrentInjectionGroupControl(const std::string& groupName, const Group2::InjectionCMode& groupControl ) {
current_injection_group_controls_[groupName] = groupControl;
}
const Group2::InjectionCMode& currentInjectionGroupControl(const std::string& groupName) const {
auto it = current_injection_group_controls_.find(groupName);
if (it == current_injection_group_controls_.end())
OPM_THROW(std::logic_error, "Could not find any well control for " << groupName);
return it->second;
}
data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override
{
@ -538,8 +392,8 @@ namespace Opm
for( const auto& wt : this->wellMap() ) {
const auto w = wt.second[ 0 ];
auto& well = res.at( wt.first );
well.control = this->currentControls()[ w ];
well.injectionControl = static_cast<int>(this->currentInjectionControls()[ w ]);
well.productionControl = static_cast<int>(this->currentProductionControls()[ w ]);
const int well_rate_index = w * pu.num_phases;
if ( pu.phase_used[Water] ) {
@ -835,6 +689,11 @@ namespace Opm
return well_reservoir_rates_;
}
const std::vector<double>& wellReservoirRates() const
{
return well_reservoir_rates_;
}
std::vector<double>& wellDissolvedGasRates()
{
return well_dissolved_gas_rates_;
@ -919,7 +778,12 @@ namespace Opm
private:
std::vector<double> perfphaserates_;
std::vector<int> current_controls_;
std::vector<Opm::Well2::InjectorCMode> current_injection_controls_;
std::vector<Well2::ProducerCMode> current_production_controls_;
std::map<std::string, Group2::ProductionCMode> current_production_group_controls_;
std::map<std::string, Group2::InjectionCMode> current_injection_group_controls_;
std::vector<double> perfRateSolvent_;
// it is the throughput of water flow through the perforations

7
tests/test_ecl_output.cc
View File

@ -222,7 +222,9 @@ void test_readWriteWells()
w1.rates = r1;
w1.bhp = 1.23;
w1.temperature = 3.45;
w1.control = 1;
w1.injectionControl = 1;
w1.productionControl = 1;
/*
* the connection keys (active indices) and well names correspond to the
@ -234,7 +236,8 @@ void test_readWriteWells()
w2.rates = r2;
w2.bhp = 2.34;
w2.temperature = 4.56;
w2.control = 2;
w1.injectionControl = 2;
w1.productionControl = 2;
w2.connections.push_back( { 188, rc3, 36.22, 19.28, 0.0, 0.0, 0.0, 0.0 } );
Opm::data::Wells wellRates;

22
tests/test_wellmodel.cpp
View File

@ -206,16 +206,6 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
BOOST_CHECK(well->wellType() == PRODUCER);
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4);
const auto& wc = well->wellControls();
const int ctrl_num = well_controls_get_num(wc);
BOOST_CHECK(ctrl_num > 0);
const auto& control = well_controls_get_current(wc);
BOOST_CHECK(control >= 0);
// GAS RATE CONTROL
const auto& distr = well_controls_iget_distr(wc, control);
BOOST_CHECK(distr[0] == 0.);
BOOST_CHECK(distr[1] == 0.);
BOOST_CHECK(distr[2] == 1.);
}
// second well, it is the injection well from the deck
@ -224,16 +214,6 @@ BOOST_AUTO_TEST_CASE(TestBehavoir) {
BOOST_CHECK_EQUAL(well->name(), "INJE1");
BOOST_CHECK(well->wellType() == INJECTOR);
BOOST_CHECK(well->numEq == 3);
BOOST_CHECK(well->numStaticWellEq== 4);
const auto& wc = well->wellControls();
const int ctrl_num = well_controls_get_num(wc);
BOOST_CHECK(ctrl_num > 0);
const auto& control = well_controls_get_current(wc);
BOOST_CHECK(control >= 0);
// WATER RATE CONTROL
const auto& distr = well_controls_iget_distr(wc, control);
BOOST_CHECK(distr[0] == 1.);
BOOST_CHECK(distr[1] == 0.);
BOOST_CHECK(distr[2] == 0.);
BOOST_CHECK(well->numStaticWellEq== 4);
}
}