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

Merge pull request #3498 from totto82/gpmaint

Browse Source 
Support for gpmaint
This commit is contained in:
Tor Harald Sandve
2021-09-07 12:38:48 +02:00
committed by GitHub
parent 0aa60a6f7c adda160098
commit 3b5a040197
15 changed files with 957 additions and 392 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@@ -309,6 +309,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/simulators/wells/PerfData.hpp
opm/simulators/wells/PerforationData.hpp
opm/simulators/wells/RateConverter.hpp
opm/simulators/wells/RegionAttributeHelpers.hpp
opm/simulators/wells/RegionAverageCalculator.hpp
opm/simulators/utils/readDeck.hpp
opm/simulators/wells/SingleWellState.hpp
opm/simulators/wells/TargetCalculator.hpp

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

@@ -60,6 +60,7 @@
#include <opm/simulators/wells/WellState.hpp>
#include <opm/simulators/wells/WGState.hpp>
#include <opm/simulators/wells/RateConverter.hpp>
#include <opm/simulators/wells/RegionAverageCalculator.hpp>
#include <opm/simulators/wells/WellInterface.hpp>
#include <opm/simulators/wells/StandardWell.hpp>
#include <opm/simulators/wells/MultisegmentWell.hpp>
@@ -132,6 +133,10 @@ namespace Opm {
using RateConverterType = RateConverter::
SurfaceToReservoirVoidage<FluidSystem, std::vector<int> >;
// For computing average pressured used by gpmaint
using AverageRegionalPressureType = RegionAverageCalculator::
AverageRegionalPressure<FluidSystem, std::vector<int> >;
BlackoilWellModel(Simulator& ebosSimulator);
void init();
@@ -317,6 +322,8 @@ namespace Opm {
bool alternative_well_rate_init_{};
std::unique_ptr<RateConverterType> rateConverter_{};
std::unique_ptr<AverageRegionalPressureType> regionalAveragePressureCalculator_{};
SimulatorReportSingle last_report_{};

28
opm/simulators/wells/BlackoilWellModelGeneric.cpp
View File

@@ -450,10 +450,15 @@ checkGroupInjectionConstraints(const Group& group,
current_rate = comm_.sum(current_rate);
const auto& controls = group.injectionControls(phase, this->summaryState_);
if (controls.surface_max_rate < current_rate) {
double target = controls.surface_max_rate;
if (group.has_gpmaint_control(phase, Group::InjectionCMode::RATE) && this->groupState().has_gpmaint_target(group.name()))
target = this->groupState().gpmaint_target(group.name());
if (target < current_rate) {
double scale = 1.0;
if (current_rate > 1e-12)
scale = controls.surface_max_rate / current_rate;
scale = target / current_rate;
return std::make_pair(Group::InjectionCMode::RATE, scale);
}
}
@@ -468,10 +473,15 @@ checkGroupInjectionConstraints(const Group& group,
current_rate = comm_.sum(current_rate);
const auto& controls = group.injectionControls(phase, this->summaryState_);
if (controls.resv_max_rate < current_rate) {
double target = controls.resv_max_rate;
if (group.has_gpmaint_control(phase, Group::InjectionCMode::RESV) && this->groupState().has_gpmaint_target(group.name()))
target = this->groupState().gpmaint_target(group.name());
if (target < current_rate) {
double scale = 1.0;
if (current_rate > 1e-12)
scale = controls.resv_max_rate / current_rate;
scale = target / current_rate;
return std::make_pair(Group::InjectionCMode::RESV, scale);
}
}
@@ -637,10 +647,14 @@ checkGroupProductionConstraints(const Group& group,
// sum over all nodes
current_rate = comm_.sum(current_rate);
if (controls.resv_target < current_rate ) {
double target = controls.resv_target;
if (group.has_gpmaint_control(Group::ProductionCMode::RESV) && this->groupState().has_gpmaint_target(group.name()))
target = this->groupState().gpmaint_target(group.name());
if ( target < current_rate ) {
double scale = 1.0;
if (current_rate > 1e-12)
scale = controls.resv_target / current_rate;
scale = target / current_rate;
return std::make_pair(Group::ProductionCMode::RESV, scale);
}
}
@@ -1521,6 +1535,8 @@ updateAndCommunicateGroupData(const int reportStepIdx,
WellGroupHelpers::updateVREPForGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState());
WellGroupHelpers::updateReservoirRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState());
WellGroupHelpers::updateSurfaceRatesInjectionGroups(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, this->groupState());
WellGroupHelpers::updateGroupProductionRates(fieldGroup, schedule(), reportStepIdx, well_state_nupcol, well_state, this->groupState());
// We use the rates from the previous time-step to reduce oscillations

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

@@ -203,10 +203,15 @@ namespace Opm {
WellGroupHelpers::setCmodeGroup(fieldGroup, schedule(), summaryState, timeStepIdx, this->wellState(), this->groupState());
// Compute reservoir volumes for RESV controls.
rateConverter_.reset(new RateConverterType (phase_usage_,
std::vector<int>(local_num_cells_, 0)));
rateConverter_.reset(new RateConverterType (phase_usage_, std::vector<int>(local_num_cells_, 0)));
rateConverter_->template defineState<ElementContext>(ebosSimulator_);
// Compute regional average pressures used by gpmaint
const auto& fp = this->eclState_.fieldProps();
const auto& fipnum = fp.get_int("FIPNUM");
regionalAveragePressureCalculator_.reset(new AverageRegionalPressureType (phase_usage_,fipnum));
regionalAveragePressureCalculator_->template defineState<ElementContext>(ebosSimulator_);
{
const auto& sched_state = this->schedule()[timeStepIdx];
// update VFP properties
@@ -458,6 +463,12 @@ namespace Opm {
// update the rate converter with current averages pressures etc in
rateConverter_->template defineState<ElementContext>(ebosSimulator_);
regionalAveragePressureCalculator_->template defineState<ElementContext>(ebosSimulator_);
const Group& fieldGroup = schedule_.getGroup("FIELD", reportStepIdx);
WellGroupHelpers::updateGpMaintTargetForGroups(fieldGroup,
schedule_, *regionalAveragePressureCalculator_, reportStepIdx, dt, this->wellState(), this->groupState());
// calculate the well potentials
try {
updateWellPotentials(reportStepIdx,
@@ -472,7 +483,6 @@ namespace Opm {
updateWellTestState(simulationTime, wellTestState_);
// check group sales limits at the end of the timestep
const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
checkGconsaleLimits(fieldGroup, this->wellState(),
ebosSimulator_.episodeIndex(), local_deferredLogger);

39
opm/simulators/wells/GroupState.cpp
View File

@@ -107,6 +107,27 @@ const std::vector<double>& GroupState::injection_reduction_rates(const std::stri
return group_iter->second;
}
//-------------------------------------------------------------------------
bool GroupState::has_injection_surface_rates(const std::string& gname) const {
auto group_iter = this->inj_surface_rates.find(gname);
return (group_iter != this->inj_surface_rates.end());
}
void GroupState::update_injection_surface_rates(const std::string& gname, const std::vector<double>& rates) {
if (rates.size() != this->num_phases)
throw std::logic_error("Wrong number of phases");
this->inj_surface_rates[gname] = rates;
}
const std::vector<double>& GroupState::injection_surface_rates(const std::string& gname) const {
auto group_iter = this->inj_surface_rates.find(gname);
if (group_iter == this->inj_surface_rates.end())
throw std::logic_error("No such group");
return group_iter->second;
}
//-------------------------------------------------------------------------
@@ -246,6 +267,24 @@ GPMaint::State& GroupState::gpmaint(const std::string& gname) {
}
//-------------------------------------------------------------------------
void GroupState::update_gpmaint_target(const std::string& gname, double target) {
this->m_gpmaint_target[gname] = target;
}
double GroupState::gpmaint_target(const std::string& gname) const {
auto group_iter = this->m_gpmaint_target.find(gname);
if (group_iter == this->m_gpmaint_target.end())
throw std::logic_error("No such group");
return group_iter->second;
}
bool GroupState::has_gpmaint_target(const std::string& gname) const {
return (this->m_gpmaint_target.count(gname) > 0);
}
//-------------------------------------------------------------------------
namespace {

11
opm/simulators/wells/GroupState.hpp
View File

@@ -52,6 +52,11 @@ public:
void update_injection_reservoir_rates(const std::string& gname, const std::vector<double>& rates);
const std::vector<double>& injection_reservoir_rates(const std::string& gname) const;
bool has_injection_surface_rates(const std::string& gname) const;
void update_injection_surface_rates(const std::string& gname, const std::vector<double>& rates);
const std::vector<double>& injection_surface_rates(const std::string& gname) const;
void update_injection_rein_rates(const std::string& gname, const std::vector<double>& rates);
const std::vector<double>& injection_rein_rates(const std::string& gname) const;
@@ -65,6 +70,10 @@ public:
double grat_sales_target(const std::string& gname) const;
bool has_grat_sales_target(const std::string& gname) const;
void update_gpmaint_target(const std::string& gname, double target);
double gpmaint_target(const std::string& gname) const;
bool has_gpmaint_target(const std::string& gname) const;
bool has_production_control(const std::string& gname) const;
void production_control(const std::string& gname, Group::ProductionCMode cmode);
Group::ProductionCMode production_control(const std::string& gname) const;
@@ -161,11 +170,13 @@ private:
std::map<std::string, Group::ProductionCMode> production_controls;
std::map<std::string, std::vector<double>> prod_red_rates;
std::map<std::string, std::vector<double>> inj_red_rates;
std::map<std::string, std::vector<double>> inj_surface_rates;
std::map<std::string, std::vector<double>> inj_resv_rates;
std::map<std::string, std::vector<double>> inj_potentials;
std::map<std::string, std::vector<double>> inj_rein_rates;
std::map<std::string, double> inj_vrep_rate;
std::map<std::string, double> m_grat_sales_target;
std::map<std::string, double> m_gpmaint_target;
std::map<std::pair<Phase, std::string>, Group::InjectionCMode> injection_controls;

403
opm/simulators/wells/RateConverter.hpp
View File

@@ -25,7 +25,7 @@
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/grid/utility/RegionMapping.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
#include <dune/grid/common/gridenums.hh>
#include <algorithm>
#include <cmath>
@@ -47,345 +47,6 @@
namespace Opm {
namespace RateConverter {
/**
* Convenience tools for implementing the rate conversion
* facility.
*/
namespace Details {
namespace Select {
template <class RegionID, bool>
struct RegionIDParameter
{
using type =
typename std::remove_reference<RegionID>::type &;
};
template <class RegionID>
struct RegionIDParameter<RegionID, true>
{
using type = RegionID;
};
} // Select
/**
* \brief Computes the temperature, pressure, and counter increment.
*
* In a parallel run only cells owned contribute to the cell average.
* \tparam is_parallel Whether this is a parallel run.
*/
template<bool is_parallel>
struct AverageIncrementCalculator
{
/**
* \brief Computes the temperature, pressure, and counter increment.
* \param pressure The pressure.
* \param temperature The temperature.
* \param rs The rs.
* \param rv The rv.
* \param cell The current cell index.
* \param ownership A vector indicating whether a cell is owned
* by this process (value 1), or not (value 0).
* \param cell The cell index.
*/
std::tuple<double, double, double, double, int>
operator()(const std::vector<double>& pressure,
const std::vector<double>& temperature,
const std::vector<double>& rs,
const std::vector<double>& rv,
const std::vector<double>& ownership,
std::size_t cell){
if ( ownership[cell] )
{
return std::make_tuple(pressure[cell],
temperature[cell],
rs[cell],
rv[cell],
1);
}
else
{
return std::make_tuple(0, 0, 0, 0, 0);
}
}
};
template<>
struct AverageIncrementCalculator<false>
{
std::tuple<double, double, double, double, int>
operator()(const std::vector<double>& pressure,
const std::vector<double>& temperature,
const std::vector<double>& rs,
const std::vector<double>& rv,
const std::vector<double>&,
std::size_t cell){
return std::make_tuple(pressure[cell],
temperature[cell],
rs[cell],
rv[cell],
1);
}
};
/**
* Provide mapping from Region IDs to user-specified collection
* of per-region attributes.
*
* \tparam RegionId Region identifier type. Must be hashable by
* \code std::hash<> \endcode. Typically a built-in integer
* type--e.g., \c int.
*
* \tparam Attributes User-defined type that represents
* collection of attributes that have meaning in a per-region
* aggregate sense. Must be copy-constructible.
*/
template <typename RegionId, class Attributes>
class RegionAttributes
{
public:
/**
* Expose \c RegionId as a vocabulary type for use in query
* methods.
*/
using RegionID =
typename Select::RegionIDParameter
<RegionId, std::is_integral<RegionId>::value>::type;
/**
* Constructor.
*
* \tparam RMap Class type that implements the RegionMapping
* protocol. Typically an instantiation of \code
* Opm::RegionMapping<> \endcode.
*
* \param[in] rmap Specific region mapping that provides
* reverse lookup from regions to cells.
*
* \param[in] attr Pre-constructed initialiser for \c
* Attributes.
*/
template <class RMap>
RegionAttributes(const RMap& rmap,
const Attributes& attr)
{
using VT = typename AttributeMap::value_type;
for (const auto& r : rmap.activeRegions()) {
auto v = std::make_unique<Value>(attr);
const auto stat = attr_.insert(VT(r, std::move(v)));
if (stat.second) {
// New value inserted.
const auto& cells = rmap.cells(r);
assert (! cells.empty());
// Region's representative cell.
stat.first->second->cell_ = cells[0];
}
}
}
/**
* Retrieve representative cell in region.
*
* \param[in] reg Specific region.
*
* \return Representative cell in region \p reg.
*/
int cell(const RegionID reg) const
{
return this->find(reg).cell_;
}
/**
* Request read-only access to region's attributes.
*
* \param[in] reg Specific region.
*
* \return Read-only access to region \p reg's per-region
* attributes.
*/
const Attributes& attributes(const RegionID reg) const
{
return this->find(reg).attr_;
}
/**
* Request modifiable access to region's attributes.
*
* \param[in] reg Specific region.
*
* \return Read-write access to region \p reg's per-region
* attributes.
*/
Attributes& attributes(const RegionID reg)
{
return this->find(reg).attr_;
}
private:
/**
* Aggregate per-region attributes along with region's
* representative cell.
*/
struct Value {
Value(const Attributes& attr)
: attr_(attr)
, cell_(-1)
{}
Attributes attr_;
int cell_;
};
using ID =
typename std::remove_reference<RegionId>::type;
using AttributeMap =
std::unordered_map<ID, std::unique_ptr<Value>>;
AttributeMap attr_;
/**
* Read-only access to region's properties.
*/
const Value& find(const RegionID reg) const
{
const auto& i = attr_.find(reg);
if (i == attr_.end()) {
throw std::invalid_argument("Unknown region ID");
}
return *i->second;
}
/**
* Read-write access to region's properties.
*/
Value& find(const RegionID reg)
{
const auto& i = attr_.find(reg);
if (i == attr_.end()) {
throw std::invalid_argument("Unknown region ID");
}
return *i->second;
}
};
/**
* Convenience functions for querying presence/absence of
* active phases.
*/
namespace PhaseUsed {
/**
* Active water predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not water is an active phase.
*/
inline bool
water(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
}
/**
* Active oil predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not oil is an active phase.
*/
inline bool
oil(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
}
/**
* Active gas predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not gas is an active phase.
*/
inline bool
gas(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
}
} // namespace PhaseUsed
/**
* Convenience functions for querying numerical IDs
* ("positions") of active phases.
*/
namespace PhasePos {
/**
* Numerical ID of active water phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of water if
* active, -1 if not.
*/
inline int
water(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::water(pu)) {
p = pu.phase_pos[ BlackoilPhases::Aqua ];
}
return p;
}
/**
* Numerical ID of active oil phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of oil if
* active, -1 if not.
*/
inline int
oil(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::oil(pu)) {
p = pu.phase_pos[ BlackoilPhases::Liquid ];
}
return p;
}
/**
* Numerical ID of active gas phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of gas if
* active, -1 if not.
*/
inline int
gas(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::gas(pu)) {
p = pu.phase_pos[ BlackoilPhases::Vapour ];
}
return p;
}
} // namespace PhasePos
} // namespace Details
/**
* Convert component rates at surface conditions to phase
@@ -484,7 +145,7 @@ namespace Opm {
// only count oil and gas filled parts of the domain
double hydrocarbon = 1.0;
const auto& pu = phaseUsage_;
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
}
@@ -496,15 +157,15 @@ namespace Opm {
if (hydrocarbonPV > 0.) {
auto& attr = attributes_hpv[reg];
attr.pv += hydrocarbonPV;
if (Details::PhaseUsed::oil(pu) && Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu) && RegionAttributeHelpers::PhaseUsed::gas(pu)) {
attr.rs += fs.Rs().value() * hydrocarbonPV;
attr.rv += fs.Rv().value() * hydrocarbonPV;
}
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
} else {
assert(Details::PhaseUsed::gas(pu));
assert(RegionAttributeHelpers::PhaseUsed::gas(pu));
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
}
@@ -514,18 +175,18 @@ namespace Opm {
if (pv_cell > 0.) {
auto& attr = attributes_pv[reg];
attr.pv += pv_cell;
if (Details::PhaseUsed::oil(pu) && Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu) && RegionAttributeHelpers::PhaseUsed::gas(pu)) {
attr.rs += fs.Rs().value() * pv_cell;
attr.rv += fs.Rv().value() * pv_cell;
}
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
attr.temperature += fs.temperature(FluidSystem::oilPhaseIdx).value() * pv_cell;
} else if (Details::PhaseUsed::gas(pu)) {
} else if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
attr.temperature += fs.temperature(FluidSystem::gasPhaseIdx).value() * pv_cell;
} else {
assert(Details::PhaseUsed::water(pu));
assert(RegionAttributeHelpers::PhaseUsed::water(pu));
attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
attr.temperature += fs.temperature(FluidSystem::waterPhaseIdx).value() * pv_cell;
}
@@ -616,13 +277,13 @@ namespace Opm {
const double T = ra.temperature;
const double saltConcentration = ra.saltConcentration;
const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu);
const int iw = RegionAttributeHelpers::PhasePos::water(pu);
const int io = RegionAttributeHelpers::PhasePos::oil (pu);
const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
@@ -637,7 +298,7 @@ namespace Opm {
// Determinant of 'R' matrix
const double detR = 1.0 - (Rs * Rv);
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
// q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
@@ -645,12 +306,12 @@ namespace Opm {
coeff[io] += 1.0 / den;
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
coeff[ig] -= ra.rv / den;
}
}
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
@@ -658,7 +319,7 @@ namespace Opm {
coeff[ig] += 1.0 / den;
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
coeff[io] -= ra.rs / den;
}
}
@@ -674,13 +335,13 @@ namespace Opm {
const double T = ra.temperature;
const double saltConcentration = ra.saltConcentration;
const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu);
const int iw = RegionAttributeHelpers::PhasePos::water(pu);
const int io = RegionAttributeHelpers::PhasePos::oil (pu);
const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
std::fill(& coeff[0], & coeff[0] + phaseUsage_.num_phases, 0.0);
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
@@ -688,12 +349,12 @@ namespace Opm {
coeff[iw] = 1.0 / bw;
}
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
coeff[io] += 1.0 / bo;
}
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, 0.0);
coeff[ig] += 1.0 / bg;
}
@@ -732,11 +393,11 @@ namespace Opm {
const double T = ra.temperature;
const double saltConcentration = ra.saltConcentration;
const int iw = Details::PhasePos::water(pu);
const int io = Details::PhasePos::oil (pu);
const int ig = Details::PhasePos::gas (pu);
const int iw = RegionAttributeHelpers::PhasePos::water(pu);
const int io = RegionAttributeHelpers::PhasePos::oil (pu);
const int ig = RegionAttributeHelpers::PhasePos::gas (pu);
if (Details::PhaseUsed::water(pu)) {
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
// q[w]_r = q[w]_s / bw
const double bw = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, saltConcentration);
@@ -764,7 +425,7 @@ namespace Opm {
// Determinant of 'R' matrix
const double detR = 1.0 - (Rs * Rv);
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
// q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)
const double bo = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
@@ -772,14 +433,14 @@ namespace Opm {
voidage_rates[io] = surface_rates[io];
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
voidage_rates[io] -= Rv * surface_rates[ig];
}
voidage_rates[io] /= den;
}
if (Details::PhaseUsed::gas(pu)) {
if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
// q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)
const double bg = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
@@ -787,7 +448,7 @@ namespace Opm {
voidage_rates[ig] = surface_rates[ig];
if (Details::PhaseUsed::oil(pu)) {
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
voidage_rates[ig] -= Rs * surface_rates[io];
}
@@ -852,7 +513,7 @@ namespace Opm {
double saltConcentration;
};
Details::RegionAttributes<RegionId, Attributes> attr_;
RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
};
} // namespace RateConverter

410
opm/simulators/wells/RegionAttributeHelpers.hpp Normal file
View File

@@ -0,0 +1,410 @@
/*
Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
Copyright 2014, 2015 Statoil ASA.
Copyright 2017, IRIS
Copyright 2017, Equinor
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_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED
#define OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/grid/utility/RegionMapping.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <dune/grid/common/gridenums.hh>
#include <algorithm>
#include <cmath>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>
namespace Opm {
namespace RegionAttributeHelpers {
/**
* Convenience tools for processing region
* spesific attributes
*/
namespace Select {
template <class RegionID, bool>
struct RegionIDParameter
{
using type =
typename std::remove_reference<RegionID>::type &;
};
template <class RegionID>
struct RegionIDParameter<RegionID, true>
{
using type = RegionID;
};
} // Select
/**
* \brief Computes the temperature, pressure, and counter increment.
*
* In a parallel run only cells owned contribute to the cell average.
* \tparam is_parallel Whether this is a parallel run.
*/
template<bool is_parallel>
struct AverageIncrementCalculator
{
/**
* \brief Computes the temperature, pressure, and counter increment.
* \param pressure The pressure.
* \param temperature The temperature.
* \param rs The rs.
* \param rv The rv.
* \param cell The current cell index.
* \param ownership A vector indicating whether a cell is owned
* by this process (value 1), or not (value 0).
* \param cell The cell index.
*/
std::tuple<double, double, double, double, int>
operator()(const std::vector<double>& pressure,
const std::vector<double>& temperature,
const std::vector<double>& rs,
const std::vector<double>& rv,
const std::vector<double>& ownership,
std::size_t cell){
if ( ownership[cell] )
{
return std::make_tuple(pressure[cell],
temperature[cell],
rs[cell],
rv[cell],
1);
}
else
{
return std::make_tuple(0, 0, 0, 0, 0);
}
}
};
template<>
struct AverageIncrementCalculator<false>
{
std::tuple<double, double, double, double, int>
operator()(const std::vector<double>& pressure,
const std::vector<double>& temperature,
const std::vector<double>& rs,
const std::vector<double>& rv,
const std::vector<double>&,
std::size_t cell){
return std::make_tuple(pressure[cell],
temperature[cell],
rs[cell],
rv[cell],
1);
}
};
/**
* Provide mapping from Region IDs to user-specified collection
* of per-region attributes.
*
* \tparam RegionId Region identifier type. Must be hashable by
* \code std::hash<> \endcode. Typically a built-in integer
* type--e.g., \c int.
*
* \tparam Attributes User-defined type that represents
* collection of attributes that have meaning in a per-region
* aggregate sense. Must be copy-constructible.
*/
template <typename RegionId, class Attributes>
class RegionAttributes
{
public:
/**
* Expose \c RegionId as a vocabulary type for use in query
* methods.
*/
using RegionID =
typename Select::RegionIDParameter
<RegionId, std::is_integral<RegionId>::value>::type;
using ID =
typename std::remove_reference<RegionId>::type;
/**
* Aggregate per-region attributes along with region's
* representative cell.
*/
struct Value {
Value(const Attributes& attr)
: attr_(attr)
, cell_(-1)
{}
Attributes attr_;
int cell_;
};
using AttributeMap =
std::unordered_map<ID, std::unique_ptr<Value>>;
/**
* Constructor.
*
* \tparam RMap Class type that implements the RegionMapping
* protocol. Typically an instantiation of \code
* Opm::RegionMapping<> \endcode.
*
* \param[in] rmap Specific region mapping that provides
* reverse lookup from regions to cells.
*
* \param[in] attr Pre-constructed initialiser for \c
* Attributes.
*/
template <class RMap>
RegionAttributes(const RMap& rmap,
const Attributes& attr)
{
using VT = typename AttributeMap::value_type;
for (const auto& r : rmap.activeRegions()) {
auto v = std::make_unique<Value>(attr);
const auto stat = attr_.insert(VT(r, std::move(v)));
if (stat.second) {
// New value inserted.
const auto& cells = rmap.cells(r);
assert (! cells.empty());
// Region's representative cell.
stat.first->second->cell_ = cells[0];
}
}
}
/**
* Retrieve representative cell in region.
*
* \param[in] reg Specific region.
*
* \return Representative cell in region \p reg.
*/
int cell(const RegionID reg) const
{
return this->find(reg).cell_;
}
bool has(const RegionID reg) const
{
return this->attr_.find(reg) != this->attr_.end();
}
void insert(const RegionID r, const Attributes& attr)
{
auto [pos, inserted] = this->attr_.try_emplace(r, std::make_unique<Value>(attr));
if (inserted) {
pos->second->cell_ = -1; // NOT -1.0 -- "cell_" is 'int'
}
}
/**
* Request read-only access to region's attributes.
*
*
* \return Read-only access to all regions attributes.
*/
const AttributeMap& attributes() const
{
return attr_;
}
/**
* Request read-only access to region's attributes.
*
* \param[in] reg Specific region.
*
* \return Read-only access to region \p reg's per-region
* attributes.
*/
const Attributes& attributes(const RegionID reg) const
{
return this->find(reg).attr_;
}
/**
* Request modifiable access to region's attributes.
*
* \param[in] reg Specific region.
*
* \return Read-write access to region \p reg's per-region
* attributes.
*/
Attributes& attributes(const RegionID reg)
{
return this->find(reg).attr_;
}
private:
AttributeMap attr_;
/**
* Read-only access to region's properties.
*/
const Value& find(const RegionID reg) const
{
const auto& i = attr_.find(reg);
if (i == attr_.end()) {
throw std::invalid_argument("Unknown region ID");
}
return *i->second;
}
/**
* Read-write access to region's properties.
*/
Value& find(const RegionID reg)
{
const auto& i = attr_.find(reg);
if (i == attr_.end()) {
throw std::invalid_argument("Unknown region ID");
}
return *i->second;
}
};
/**
* Convenience functions for querying presence/absence of
* active phases.
*/
namespace PhaseUsed {
/**
* Active water predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not water is an active phase.
*/
inline bool
water(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
}
/**
* Active oil predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not oil is an active phase.
*/
inline bool
oil(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
}
/**
* Active gas predicate.
*
* \param[in] pu Active phase object.
*
* \return Whether or not gas is an active phase.
*/
inline bool
gas(const PhaseUsage& pu)
{
return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
}
} // namespace PhaseUsed
/**
* Convenience functions for querying numerical IDs
* ("positions") of active phases.
*/
namespace PhasePos {
/**
* Numerical ID of active water phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of water if
* active, -1 if not.
*/
inline int
water(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::water(pu)) {
p = pu.phase_pos[ BlackoilPhases::Aqua ];
}
return p;
}
/**
* Numerical ID of active oil phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of oil if
* active, -1 if not.
*/
inline int
oil(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::oil(pu)) {
p = pu.phase_pos[ BlackoilPhases::Liquid ];
}
return p;
}
/**
* Numerical ID of active gas phase.
*
* \param[in] pu Active phase object.
*
* \return Non-negative index/position of gas if
* active, -1 if not.
*/
inline int
gas(const PhaseUsage& pu)
{
int p = -1;
if (PhaseUsed::gas(pu)) {
p = pu.phase_pos[ BlackoilPhases::Vapour ];
}
return p;
}
} // namespace PhasePos
} // namespace RegionAttributesHelpers
} // namespace Opm
#endif /* OPM_REGIONATTRIBUTEHELPERS_HPP_HEADER_INCLUDED */

260
opm/simulators/wells/RegionAverageCalculator.hpp Normal file
View File

@@ -0,0 +1,260 @@
/*
Copyright 2021, Equinor
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_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
#define OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/simulators/wells/RegionAttributeHelpers.hpp>
#include <opm/simulators/linalg/ParallelIstlInformation.hpp>
#include <dune/grid/common/gridenums.hh>
#include <algorithm>
#include <cmath>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>
/**
* \file
* Facility for converting component rates at surface conditions to
* phase (voidage) rates at reservoir conditions.
*
* This uses the average hydrocarbon pressure to define fluid
* properties. The facility is intended to support Reservoir Voidage
* rates only ('RESV').
*/
namespace Opm {
namespace RegionAverageCalculator {
/**
* Computes hydrocarbon weighed average pressures over regions
*
* \tparam FluidSystem Fluid system class. Expected to be a BlackOilFluidSystem
*
* \tparam Region Type of a forward region mapping. Expected
* to provide indexed access through \code operator[]()
* \endcode as well as inner types \c value_type, \c
* size_type, and \c const_iterator. Typically \code
* std::vector<int> \endcode.
*/
template <class FluidSystem, class Region>
class AverageRegionalPressure {
public:
/**
* Constructor.
*
* \param[in] region Forward region mapping. Often
* corresponds to the "FIPNUM" mapping of an ECLIPSE input
* deck.
*/
AverageRegionalPressure(const PhaseUsage& phaseUsage,
const Region& region)
: phaseUsage_(phaseUsage)
, rmap_ (region)
, attr_ (rmap_, Attributes())
{
}
/**
* Compute pore volume averaged hydrocarbon state pressure, *
*/
template <typename ElementContext, class EbosSimulator>
void defineState(const EbosSimulator& simulator)
{
int numRegions = 0;
const auto& gridView = simulator.gridView();
const auto& comm = gridView.comm();
for (const auto& reg : rmap_.activeRegions()) {
numRegions = std::max(numRegions, reg);
}
numRegions = comm.max(numRegions);
for (int reg = 1; reg <= numRegions ; ++ reg) {
if(!attr_.has(reg))
attr_.insert(reg, Attributes());
}
// create map from cell to region
// and set all attributes to zero
for (int reg = 1; reg <= numRegions ; ++ reg) {
auto& ra = attr_.attributes(reg);
ra.pressure = 0.0;
ra.pv = 0.0;
}
// quantities for pore volume average
std::unordered_map<RegionId, Attributes> attributes_pv;
// quantities for hydrocarbon volume average
std::unordered_map<RegionId, Attributes> attributes_hpv;
for (int reg = 1; reg <= numRegions ; ++ reg) {
attributes_pv.insert({reg, Attributes()});
attributes_hpv.insert({reg, Attributes()});
}
ElementContext elemCtx( simulator );
const auto& elemEndIt = gridView.template end</*codim=*/0>();
for (auto elemIt = gridView.template begin</*codim=*/0>();
elemIt != elemEndIt;
++elemIt)
{
const auto& elem = *elemIt;
if (elem.partitionType() != Dune::InteriorEntity)
continue;
elemCtx.updatePrimaryStencil(elem);
elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
const unsigned cellIdx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& intQuants = elemCtx.intensiveQuantities(/*spaceIdx=*/0, /*timeIdx=*/0);
const auto& fs = intQuants.fluidState();
// use pore volume weighted averages.
const double pv_cell =
simulator.model().dofTotalVolume(cellIdx)
* intQuants.porosity().value();
// only count oil and gas filled parts of the domain
double hydrocarbon = 1.0;
const auto& pu = phaseUsage_;
if (RegionAttributeHelpers::PhaseUsed::water(pu)) {
hydrocarbon -= fs.saturation(FluidSystem::waterPhaseIdx).value();
}
const int reg = rmap_.region(cellIdx);
assert(reg >= 0);
// sum p, rs, rv, and T.
const double hydrocarbonPV = pv_cell*hydrocarbon;
if (hydrocarbonPV > 0.) {
auto& attr = attributes_hpv[reg];
attr.pv += hydrocarbonPV;
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * hydrocarbonPV;
} else {
assert(RegionAttributeHelpers::PhaseUsed::gas(pu));
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * hydrocarbonPV;
}
}
if (pv_cell > 0.) {
auto& attr = attributes_pv[reg];
attr.pv += pv_cell;
if (RegionAttributeHelpers::PhaseUsed::oil(pu)) {
attr.pressure += fs.pressure(FluidSystem::oilPhaseIdx).value() * pv_cell;
} else if (RegionAttributeHelpers::PhaseUsed::gas(pu)) {
attr.pressure += fs.pressure(FluidSystem::gasPhaseIdx).value() * pv_cell;
} else {
assert(RegionAttributeHelpers::PhaseUsed::water(pu));
attr.pressure += fs.pressure(FluidSystem::waterPhaseIdx).value() * pv_cell;
}
}
}
for (int reg = 1; reg <= numRegions ; ++ reg) {
auto& ra = attr_.attributes(reg);
const double hpv_sum = comm.sum(attributes_hpv[reg].pv);
// TODO: should we have some epsilon here instead of zero?
if (hpv_sum > 0.) {
const auto& attri_hpv = attributes_hpv[reg];
const double p_hpv_sum = comm.sum(attri_hpv.pressure);
ra.pressure = p_hpv_sum / hpv_sum;
} else {
// using the pore volume to do the averaging
const auto& attri_pv = attributes_pv[reg];
const double pv_sum = comm.sum(attri_pv.pv);
assert(pv_sum > 0.);
const double p_pv_sum = comm.sum(attri_pv.pressure);
ra.pressure = p_pv_sum / pv_sum;
}
}
}
/**
* Region identifier.
*
* Integral type.
*/
typedef typename RegionMapping<Region>::RegionId RegionId;
/**
* Average pressure
*
*/
double
pressure(const RegionId r) const
{
if (r == 0 ) // region 0 is the whole field
{
double pressure = 0.0;
int num_active_regions = 0;
for (const auto& attr : attr_.attributes()) {
const auto& value = *attr.second;
const auto& ra = value.attr_;
pressure += ra.pressure;
num_active_regions ++;
}
return pressure / num_active_regions;
}
const auto& ra = attr_.attributes(r);
return ra.pressure;
}
private:
/**
* Fluid property object.
*/
const PhaseUsage phaseUsage_;
/**
* "Fluid-in-place" region mapping (forward and reverse).
*/
const RegionMapping<Region> rmap_;
/**
* Derived property attributes for each active region.
*/
struct Attributes {
Attributes()
: pressure(0.0)
, pv(0.0)
{}
double pressure;
double pv;
};
RegionAttributeHelpers::RegionAttributes<RegionId, Attributes> attr_;
};
} // namespace RegionAverageCalculator
} // namespace Opm
#endif /* OPM_REGIONAVERAGECALCULATOR_HPP_HEADER_INCLUDED */

24
opm/simulators/wells/TargetCalculator.cpp
View File

@@ -40,13 +40,19 @@ namespace WellGroupHelpers
TargetCalculator::TargetCalculator(const Group::ProductionCMode cmode,
const PhaseUsage& pu,
const std::vector<double>& resv_coeff,
const double group_grat_target_from_sales)
const double group_grat_target_from_sales,
const std::string& group_name,
const GroupState& group_state,
const bool use_gpmaint)
: cmode_(cmode)
, pu_(pu)
, resv_coeff_(resv_coeff)
, group_grat_target_from_sales_(group_grat_target_from_sales)
{
, group_name_(group_name)
, group_state_(group_state)
, use_gpmaint_(use_gpmaint)
{
}
template <typename RateType>
@@ -107,7 +113,12 @@ double TargetCalculator::groupTarget(const Group::ProductionControls ctrl) const
case Group::ProductionCMode::LRAT:
return ctrl.liquid_target;
case Group::ProductionCMode::RESV:
{
if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_))
return this->group_state_.gpmaint_target(this->group_name_);
return ctrl.resv_target;
}
default:
// Should never be here.
assert(false);
@@ -142,6 +153,7 @@ InjectionTargetCalculator::InjectionTargetCalculator(const Group::InjectionCMode
const double sales_target,
const GroupState& group_state,
const Phase& injection_phase,
const bool use_gpmaint,
DeferredLogger& deferred_logger)
: cmode_(cmode)
, pu_(pu)
@@ -149,6 +161,8 @@ InjectionTargetCalculator::InjectionTargetCalculator(const Group::InjectionCMode
, group_name_(group_name)
, sales_target_(sales_target)
, group_state_(group_state)
, use_gpmaint_(use_gpmaint)
{
// initialize to avoid warning
pos_ = pu.phase_pos[BlackoilPhases::Aqua];
@@ -182,8 +196,14 @@ double InjectionTargetCalculator::groupTarget(const Group::InjectionControls& ct
{
switch (cmode_) {
case Group::InjectionCMode::RATE:
if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_))
return this->group_state_.gpmaint_target(this->group_name_);
return ctrl.surface_max_rate;
case Group::InjectionCMode::RESV:
if(use_gpmaint_ && this->group_state_.has_gpmaint_target(this->group_name_))
return this->group_state_.gpmaint_target(this->group_name_) / resv_coeff_[pos_];
return ctrl.resv_max_rate / resv_coeff_[pos_];
case Group::InjectionCMode::REIN: {
double production_rate = this->group_state_.injection_rein_rates(ctrl.reinj_group)[pos_];

10
opm/simulators/wells/TargetCalculator.hpp
View File

@@ -45,7 +45,10 @@ namespace WellGroupHelpers
TargetCalculator(const Group::ProductionCMode cmode,
const PhaseUsage& pu,
const std::vector<double>& resv_coeff,
const double group_grat_target_from_sales);
const double group_grat_target_from_sales,
const std::string& group_name,
const GroupState& group_state,
const bool use_gpmaint);
template <typename RateType>
RateType calcModeRateFromRates(const std::vector<RateType>& rates) const
@@ -65,6 +68,9 @@ namespace WellGroupHelpers
const PhaseUsage& pu_;
const std::vector<double>& resv_coeff_;
const double group_grat_target_from_sales_;
const std::string& group_name_;
const GroupState& group_state_;
bool use_gpmaint_;
};
/// Based on a group control mode, extract or calculate rates, and
@@ -79,6 +85,7 @@ namespace WellGroupHelpers
const double sales_target,
const GroupState& group_state,
const Phase& injection_phase,
const bool use_gpmaint,
DeferredLogger& deferred_logger);
template <typename RateVec>
@@ -98,6 +105,7 @@ namespace WellGroupHelpers
const std::string& group_name_;
double sales_target_;
const GroupState& group_state_;
bool use_gpmaint_;
int pos_;
GuideRateModel::Target target_;
};

42
opm/simulators/wells/WellGroupHelpers.cpp
View File

@@ -21,13 +21,10 @@
#include <config.h>
#include <opm/simulators/wells/WellGroupHelpers.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group/GConSump.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group/GConSale.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group/Group.hpp>
#include <opm/simulators/utils/DeferredLogger.hpp>
#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
#include <opm/simulators/wells/GroupState.hpp>
#include <opm/simulators/wells/TargetCalculator.hpp>
#include <opm/simulators/wells/VFPProdProperties.hpp>
#include <opm/simulators/wells/WellState.hpp>
@@ -163,6 +160,17 @@ namespace WellGroupHelpers
group_state.production_control(group.name(), controls.cmode);
}
if (group.has_gpmaint_control(Group::ProductionCMode::RESV)) {
group_state.production_control(group.name(), Group::ProductionCMode::RESV);
}
for (Phase phase : all) {
if (group.has_gpmaint_control(phase, Group::InjectionCMode::RATE)) {
group_state.injection_control(group.name(), phase, Group::InjectionCMode::RATE);
} else if (group.has_gpmaint_control(phase, Group::InjectionCMode::RESV)) {
group_state.injection_control(group.name(), phase, Group::InjectionCMode::RESV);
}
}
if (schedule[reportStepIdx].gconsale().has(group.name())) {
group_state.injection_control(group.name(), Phase::GAS, Group::InjectionCMode::SALE);
}
@@ -546,6 +554,30 @@ namespace WellGroupHelpers
group_state.update_injection_reservoir_rates(group.name(), resv);
}
void updateSurfaceRatesInjectionGroups(const Group& group,
const Schedule& schedule,
const int reportStepIdx,
const WellState& wellStateNupcol,
GroupState& group_state)
{
for (const std::string& groupName : group.groups()) {
const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
updateSurfaceRatesInjectionGroups(groupTmp, schedule, reportStepIdx, wellStateNupcol, group_state);
}
const int np = wellStateNupcol.numPhases();
std::vector<double> rates(np, 0.0);
for (int phase = 0; phase < np; ++phase) {
rates[phase] = sumWellPhaseRates(false,
group,
schedule,
wellStateNupcol,
reportStepIdx,
phase,
/*isInjector*/ true);
}
group_state.update_injection_surface_rates(group.name(), rates);
}
void updateWellRates(const Group& group,
const Schedule& schedule,
const int reportStepIdx,
@@ -1118,7 +1150,7 @@ namespace WellGroupHelpers
if (group_state.has_grat_sales_target(group.name()))
gratTargetFromSales = group_state.grat_sales_target(group.name());
TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales);
TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl));
FractionCalculator fcalc(schedule, wellState, group_state, reportStepIdx, guideRate, tcalc.guideTargetMode(), pu, true, Phase::OIL);
auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); };
@@ -1245,7 +1277,7 @@ namespace WellGroupHelpers
const auto& gconsale = schedule[reportStepIdx].gconsale().get(group.name(), summaryState);
sales_target = gconsale.sales_target;
}
InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger);
InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger);
FractionCalculator fcalc(schedule, wellState, group_state, reportStepIdx, guideRate, tcalc.guideTargetMode(), pu, false, injectionPhase);
auto localFraction = [&](const std::string& child) { return fcalc.localFraction(child, name); };

89
opm/simulators/wells/WellGroupHelpers.hpp
View File

@@ -22,6 +22,11 @@
#define OPM_WELLGROUPHELPERS_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/Schedule/Group/GuideRate.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group/GPMaint.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/simulators/wells/GroupState.hpp>
#include <opm/simulators/wells/WellState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group/Group.hpp>
#include <map>
#include <string>
@@ -151,6 +156,12 @@ namespace WellGroupHelpers
WellState& wellState,
GroupState& group_state);
void updateSurfaceRatesInjectionGroups(const Group& group,
const Schedule& schedule,
const int reportStepIdx,
const WellState& wellStateNupcol,
GroupState& group_state);
void updateWellRates(const Group& group,
const Schedule& schedule,
const int reportStepIdx,
@@ -181,6 +192,84 @@ namespace WellGroupHelpers
WellState& wellState,
GroupState& group_state);
template <class RegionalValues>
void updateGpMaintTargetForGroups(const Group& group,
const Schedule& schedule,
const RegionalValues& regional_values,
const int reportStepIdx,
const double dt,
const WellState& well_state,
GroupState& group_state)
{
for (const std::string& groupName : group.groups()) {
const Group& groupTmp = schedule.getGroup(groupName, reportStepIdx);
updateGpMaintTargetForGroups(groupTmp, schedule, regional_values, reportStepIdx, dt, well_state, group_state);
}
const auto& gpm = group.gpmaint();
if(!gpm)
return;
const auto [name, number] = *gpm->region();
const double error = gpm->pressure_target() - regional_values.pressure(number);
double current_rate = 0.0;
const auto& pu = well_state.phaseUsage();
double sign = 1.0;
switch (gpm->flow_target()) {
case Opm::GPMaint::FlowTarget::RESV_PROD:
{
current_rate = -group_state.injection_vrep_rate(group.name());
sign = -1.0;
break;
}
case Opm::GPMaint::FlowTarget::RESV_OINJ:
{
if (pu.phase_used[BlackoilPhases::Liquid])
current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Liquid]];
break;
}
case Opm::GPMaint::FlowTarget::RESV_WINJ:
{
if (pu.phase_used[BlackoilPhases::Aqua])
current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Aqua]];
break;
}
case Opm::GPMaint::FlowTarget::RESV_GINJ:
{
if (pu.phase_used[BlackoilPhases::Vapour])
current_rate = group_state.injection_reservoir_rates(group.name())[pu.phase_pos[BlackoilPhases::Vapour]];
break;
}
case Opm::GPMaint::FlowTarget::SURF_OINJ:
{
if (pu.phase_used[BlackoilPhases::Liquid])
current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Liquid]];
break;
}
case Opm::GPMaint::FlowTarget::SURF_WINJ:
{
if (pu.phase_used[BlackoilPhases::Aqua])
current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Aqua]];
break;
}
case Opm::GPMaint::FlowTarget::SURF_GINJ:
{
if (pu.phase_used[BlackoilPhases::Vapour])
current_rate = group_state.injection_surface_rates(group.name())[pu.phase_pos[BlackoilPhases::Vapour]];
break;
}
default:
throw std::invalid_argument("Invalid Flow target type in GPMAINT");
}
auto& gpmaint_state = group_state.gpmaint(group.name());
double rate = gpm->rate(gpmaint_state, current_rate, error, dt);
group_state.update_gpmaint_target(group.name(), std::max(0.0, sign * rate));
}
std::map<std::string, double>
computeNetworkPressures(const Opm::Network::ExtNetwork& network,
const WellState& well_state,

4
opm/simulators/wells/WellInterfaceEval.cpp
View File

@@ -136,7 +136,7 @@ getGroupInjectionControl(const Group& group,
const auto& gconsale = schedule[baseif_.currentStep()].gconsale().get(group.name(), summaryState);
sales_target = gconsale.sales_target;
}
WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger);
WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger);
WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, false, injectionPhase);
auto localFraction = [&](const std::string& child) {
@@ -230,7 +230,7 @@ getGroupProductionControl(const Group& group,
if (group_state.has_grat_sales_target(group.name()))
gratTargetFromSales = group_state.grat_sales_target(group.name());
WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales);
WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl));
WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, baseif_.currentStep(), baseif_.guideRate(), tcalc.guideTargetMode(), pu, true, Phase::OIL);
auto localFraction = [&](const std::string& child) {

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

@@ -1002,7 +1002,7 @@ getGroupInjectionTargetRate(const Group& group,
const auto& gconsale = schedule[currentStep()].gconsale().get(group.name(), summaryState);
sales_target = gconsale.sales_target;
}
WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, deferred_logger);
WellGroupHelpers::InjectionTargetCalculator tcalc(currentGroupControl, pu, resv_coeff, group.name(), sales_target, group_state, injectionPhase, group.has_gpmaint_control(injectionPhase, currentGroupControl), deferred_logger);
WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, currentStep(), guideRate(), tcalc.guideTargetMode(), pu, false, injectionPhase);
auto localFraction = [&](const std::string& child) {
@@ -1072,7 +1072,7 @@ getGroupProductionTargetRate(const Group& group,
if (group_state.has_grat_sales_target(group.name()))
gratTargetFromSales = group_state.grat_sales_target(group.name());
WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales);
WellGroupHelpers::TargetCalculator tcalc(currentGroupControl, pu, resv_coeff, gratTargetFromSales, group.name(), group_state, group.has_gpmaint_control(currentGroupControl));
WellGroupHelpers::FractionCalculator fcalc(schedule, well_state, group_state, currentStep(), guideRate(), tcalc.guideTargetMode(), pu, true, Phase::OIL);
auto localFraction = [&](const std::string& child) {