mirror of
https://github.com/OPM/opm-simulators.git
synced 2025-02-25 18:55:30 -06:00
finishing separating the StandardWellsDense.hpp implementations.
This commit is contained in:
Kai Bao
parent
498f40f896
commit
8b38b7b8a3
@ -147,7 +147,7 @@ enum WellVariablePositions {
|
||||
|
||||
int numCells() const;
|
||||
|
||||
void resetWellControlFromState(WellState xw);
|
||||
void resetWellControlFromState(WellState xw) const;
|
||||
|
||||
const Wells& wells() const;
|
||||
|
||||
@ -243,125 +243,18 @@ enum WellVariablePositions {
|
||||
computeWellPotentials(const Simulator& ebosSimulator,
|
||||
WellState& well_state) const;
|
||||
|
||||
|
||||
WellCollection* wellCollection() const
|
||||
{
|
||||
return well_collection_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
WellCollection* wellCollection() const;
|
||||
|
||||
const std::vector<double>&
|
||||
wellPerfEfficiencyFactors() const
|
||||
{
|
||||
return well_perforation_efficiency_factors_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void calculateEfficiencyFactors()
|
||||
{
|
||||
if ( !localWellsActive() ) {
|
||||
return;
|
||||
}
|
||||
|
||||
const int nw = wells().number_of_wells;
|
||||
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
const std::string well_name = wells().name[w];
|
||||
const WellNode& well_node = wellCollection()->findWellNode(well_name);
|
||||
|
||||
const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor();
|
||||
|
||||
// assign the efficiency factor to each perforation related.
|
||||
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) {
|
||||
well_perforation_efficiency_factors_[perf] = well_efficiency_factor;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
wellPerfEfficiencyFactors() const;
|
||||
|
||||
void calculateEfficiencyFactors();
|
||||
|
||||
void computeWellVoidageRates(const WellState& well_state,
|
||||
std::vector<double>& well_voidage_rates,
|
||||
std::vector<double>& voidage_conversion_coeffs) const
|
||||
{
|
||||
if ( !localWellsActive() ) {
|
||||
return;
|
||||
}
|
||||
// TODO: for now, we store the voidage rates for all the production wells.
|
||||
// For injection wells, the rates are stored as zero.
|
||||
// We only store the conversion coefficients for all the injection wells.
|
||||
// Later, more delicate model will be implemented here.
|
||||
// And for the moment, group control can only work for serial running.
|
||||
const int nw = well_state.numWells();
|
||||
const int np = well_state.numPhases();
|
||||
|
||||
// we calculate the voidage rate for each well, that means the sum of all the phases.
|
||||
well_voidage_rates.resize(nw, 0);
|
||||
// store the conversion coefficients, while only for the use of injection wells.
|
||||
voidage_conversion_coeffs.resize(nw * np, 1.0);
|
||||
|
||||
std::vector<double> well_rates(np, 0.0);
|
||||
std::vector<double> convert_coeff(np, 1.0);
|
||||
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
const bool is_producer = wells().type[w] == PRODUCER;
|
||||
|
||||
// not sure necessary to change all the value to be positive
|
||||
if (is_producer) {
|
||||
std::transform(well_state.wellRates().begin() + np * w,
|
||||
well_state.wellRates().begin() + np * (w + 1),
|
||||
well_rates.begin(), std::negate<double>());
|
||||
|
||||
// the average hydrocarbon conditions of the whole field will be used
|
||||
const int fipreg = 0; // Not considering FIP for the moment.
|
||||
|
||||
rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff);
|
||||
well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(),
|
||||
convert_coeff.begin(), 0.0);
|
||||
} else {
|
||||
// TODO: Not sure whether will encounter situation with all zero rates
|
||||
// and whether it will cause problem here.
|
||||
std::copy(well_state.wellRates().begin() + np * w,
|
||||
well_state.wellRates().begin() + np * (w + 1),
|
||||
well_rates.begin());
|
||||
// the average hydrocarbon conditions of the whole field will be used
|
||||
const int fipreg = 0; // Not considering FIP for the moment.
|
||||
rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff);
|
||||
std::copy(convert_coeff.begin(), convert_coeff.end(),
|
||||
voidage_conversion_coeffs.begin() + np * w);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
void applyVREPGroupControl(WellState& well_state) const
|
||||
{
|
||||
if ( wellCollection()->havingVREPGroups() ) {
|
||||
std::vector<double> well_voidage_rates;
|
||||
std::vector<double> voidage_conversion_coeffs;
|
||||
computeWellVoidageRates(well_state, well_voidage_rates, voidage_conversion_coeffs);
|
||||
wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
|
||||
|
||||
// for the wells under group control, update the currentControls for the well_state
|
||||
for (const WellNode* well_node : wellCollection()->getLeafNodes()) {
|
||||
if (well_node->isInjector() && !well_node->individualControl()) {
|
||||
const int well_index = well_node->selfIndex();
|
||||
well_state.currentControls()[well_index] = well_node->groupControlIndex();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
std::vector<double>& voidage_conversion_coeffs) const;
|
||||
|
||||
void applyVREPGroupControl(WellState& well_state) const;
|
||||
|
||||
|
||||
protected:
|
||||
@ -409,210 +302,17 @@ enum WellVariablePositions {
|
||||
double dWellFractionMax() const {return param_.dwell_fraction_max_; }
|
||||
|
||||
// protected methods
|
||||
EvalWell getBhp(const int wellIdx) const {
|
||||
const WellControls* wc = wells().ctrls[wellIdx];
|
||||
if (well_controls_get_current_type(wc) == BHP) {
|
||||
EvalWell bhp = 0.0;
|
||||
const double target_rate = well_controls_get_current_target(wc);
|
||||
bhp.setValue(target_rate);
|
||||
return bhp;
|
||||
} else if (well_controls_get_current_type(wc) == THP) {
|
||||
const int control = well_controls_get_current(wc);
|
||||
const double thp = well_controls_get_current_target(wc);
|
||||
const double alq = well_controls_iget_alq(wc, control);
|
||||
const int table_id = well_controls_iget_vfp(wc, control);
|
||||
EvalWell aqua = 0.0;
|
||||
EvalWell liquid = 0.0;
|
||||
EvalWell vapour = 0.0;
|
||||
EvalWell bhp = 0.0;
|
||||
double vfp_ref_depth = 0.0;
|
||||
EvalWell getBhp(const int wellIdx) const;
|
||||
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
EvalWell getQs(const int wellIdx, const int phaseIdx) const;
|
||||
|
||||
if (active_[ Water ]) {
|
||||
aqua = getQs(wellIdx, pu.phase_pos[ Water]);
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
liquid = getQs(wellIdx, pu.phase_pos[ Oil ]);
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
vapour = getQs(wellIdx, pu.phase_pos[ Gas ]);
|
||||
}
|
||||
if (wells().type[wellIdx] == INJECTOR) {
|
||||
bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp);
|
||||
vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
|
||||
} else {
|
||||
bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq);
|
||||
vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
|
||||
}
|
||||
EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const;
|
||||
|
||||
// pick the density in the top layer
|
||||
const int perf = wells().well_connpos[wellIdx];
|
||||
const double rho = well_perforation_densities_[perf];
|
||||
const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_);
|
||||
bhp -= dp;
|
||||
return bhp;
|
||||
EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const;
|
||||
|
||||
}
|
||||
const int nw = wells().number_of_wells;
|
||||
return wellVariables_[nw*XvarWell + wellIdx];
|
||||
}
|
||||
|
||||
EvalWell getQs(const int wellIdx, const int phaseIdx) const {
|
||||
EvalWell qs = 0.0;
|
||||
const WellControls* wc = wells().ctrls[wellIdx];
|
||||
const int np = wells().number_of_phases;
|
||||
const int nw = wells().number_of_wells;
|
||||
const double target_rate = well_controls_get_current_target(wc);
|
||||
|
||||
if (wells().type[wellIdx] == INJECTOR) {
|
||||
const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx];
|
||||
if (comp_frac == 0.0)
|
||||
return qs;
|
||||
|
||||
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
|
||||
return wellVariables_[nw*XvarWell + wellIdx];
|
||||
}
|
||||
qs.setValue(target_rate);
|
||||
return qs;
|
||||
}
|
||||
|
||||
// Producers
|
||||
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) {
|
||||
return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx);
|
||||
}
|
||||
|
||||
if (well_controls_get_current_type(wc) == SURFACE_RATE) {
|
||||
// checking how many phases are included in the rate control
|
||||
// to decide wheter it is a single phase rate control or not
|
||||
const double* distr = well_controls_get_current_distr(wc);
|
||||
int num_phases_under_rate_control = 0;
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
if (distr[phase] > 0.0) {
|
||||
num_phases_under_rate_control += 1;
|
||||
}
|
||||
}
|
||||
|
||||
// there should be at least one phase involved
|
||||
assert(num_phases_under_rate_control > 0);
|
||||
|
||||
// when it is a single phase rate limit
|
||||
if (num_phases_under_rate_control == 1) {
|
||||
if (distr[phaseIdx] == 1.0) {
|
||||
qs.setValue(target_rate);
|
||||
return qs;
|
||||
}
|
||||
|
||||
int currentControlIdx = 0;
|
||||
for (int i = 0; i < np; ++i) {
|
||||
currentControlIdx += wells().comp_frac[np*wellIdx + i] * i;
|
||||
}
|
||||
|
||||
const double eps = 1e-6;
|
||||
if (wellVolumeFractionScaled(wellIdx,currentControlIdx) < eps) {
|
||||
return qs;
|
||||
}
|
||||
return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx,currentControlIdx));
|
||||
}
|
||||
|
||||
// when it is a combined two phase rate limit, such like LRAT
|
||||
// we neec to calculate the rate for the certain phase
|
||||
if (num_phases_under_rate_control == 2) {
|
||||
EvalWell combined_volume_fraction = 0.;
|
||||
for (int p = 0; p < np; ++p) {
|
||||
if (distr[p] == 1.0) {
|
||||
combined_volume_fraction += wellVolumeFractionScaled(wellIdx, p);
|
||||
}
|
||||
}
|
||||
return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / combined_volume_fraction);
|
||||
}
|
||||
|
||||
// suppose three phase combined limit is the same with RESV
|
||||
// not tested yet.
|
||||
}
|
||||
// ReservoirRate
|
||||
return target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx);
|
||||
}
|
||||
|
||||
EvalWell wellVolumeFraction(const int wellIdx, const int phaseIdx) const {
|
||||
const int nw = wells().number_of_wells;
|
||||
if (phaseIdx == Water) {
|
||||
return wellVariables_[WFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
if (phaseIdx == Gas) {
|
||||
return wellVariables_[GFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
// Oil fraction
|
||||
EvalWell well_fraction = 1.0;
|
||||
if (active_[Water]) {
|
||||
well_fraction -= wellVariables_[WFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
if (active_[Gas]) {
|
||||
well_fraction -= wellVariables_[GFrac * nw + wellIdx];
|
||||
}
|
||||
return well_fraction;
|
||||
}
|
||||
|
||||
EvalWell wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const {
|
||||
const WellControls* wc = wells().ctrls[wellIdx];
|
||||
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
|
||||
const double* distr = well_controls_get_current_distr(wc);
|
||||
return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx];
|
||||
}
|
||||
std::vector<double> g = {1,1,0.01};
|
||||
return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
template <class WellState>
|
||||
bool checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
|
||||
const WellState& well_state,
|
||||
const int well_number) const
|
||||
{
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
const int np = well_state.numPhases();
|
||||
|
||||
if (econ_production_limits.onMinOilRate()) {
|
||||
assert(active_[Oil]);
|
||||
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
|
||||
const double min_oil_rate = econ_production_limits.minOilRate();
|
||||
if (std::abs(oil_rate) < min_oil_rate) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMinGasRate() ) {
|
||||
assert(active_[Gas]);
|
||||
const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ];
|
||||
const double min_gas_rate = econ_production_limits.minGasRate();
|
||||
if (std::abs(gas_rate) < min_gas_rate) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMinLiquidRate() ) {
|
||||
assert(active_[Oil]);
|
||||
assert(active_[Water]);
|
||||
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
|
||||
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
|
||||
const double liquid_rate = oil_rate + water_rate;
|
||||
const double min_liquid_rate = econ_production_limits.minLiquidRate();
|
||||
if (std::abs(liquid_rate) < min_liquid_rate) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMinReservoirFluidRate()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
const int well_number) const;
|
||||
|
||||
using WellMapType = typename WellState::WellMapType;
|
||||
using WellMapEntryType = typename WellState::mapentry_t;
|
||||
@ -631,297 +331,18 @@ enum WellVariablePositions {
|
||||
};
|
||||
|
||||
|
||||
template <class WellState>
|
||||
RatioCheckTuple checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
|
||||
const WellState& well_state,
|
||||
const WellMapEntryType& map_entry) const
|
||||
{
|
||||
// TODO: not sure how to define the worst-offending connection when more than one
|
||||
// ratio related limit is violated.
|
||||
// The defintion used here is that we define the violation extent based on the
|
||||
// ratio between the value and the corresponding limit.
|
||||
// For each violated limit, we decide the worst-offending connection separately.
|
||||
// Among the worst-offending connections, we use the one has the biggest violation
|
||||
// extent.
|
||||
const WellMapEntryType& map_entry) const;
|
||||
|
||||
bool any_limit_violated = false;
|
||||
bool last_connection = false;
|
||||
int worst_offending_connection = INVALIDCONNECTION;
|
||||
double violation_extent = -1.0;
|
||||
|
||||
if (econ_production_limits.onMaxWaterCut()) {
|
||||
const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry);
|
||||
bool water_cut_violated = std::get<0>(water_cut_return);
|
||||
if (water_cut_violated) {
|
||||
any_limit_violated = true;
|
||||
const double violation_extent_water_cut = std::get<3>(water_cut_return);
|
||||
if (violation_extent_water_cut > violation_extent) {
|
||||
violation_extent = violation_extent_water_cut;
|
||||
worst_offending_connection = std::get<2>(water_cut_return);
|
||||
last_connection = std::get<1>(water_cut_return);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMaxGasOilRatio()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMaxWaterGasRatio()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMaxGasLiquidRatio()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
|
||||
}
|
||||
|
||||
if (any_limit_violated) {
|
||||
assert(worst_offending_connection >=0);
|
||||
assert(violation_extent > 1.);
|
||||
}
|
||||
|
||||
return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
|
||||
}
|
||||
|
||||
template <class WellState>
|
||||
RatioCheckTuple checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
|
||||
const WellState& well_state,
|
||||
const WellMapEntryType& map_entry) const
|
||||
{
|
||||
bool water_cut_limit_violated = false;
|
||||
int worst_offending_connection = INVALIDCONNECTION;
|
||||
bool last_connection = false;
|
||||
double violation_extent = -1.0;
|
||||
const WellMapEntryType& map_entry) const;
|
||||
|
||||
const int np = well_state.numPhases();
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
const int well_number = map_entry[0];
|
||||
|
||||
assert(active_[Oil]);
|
||||
assert(active_[Water]);
|
||||
|
||||
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
|
||||
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
|
||||
const double liquid_rate = oil_rate + water_rate;
|
||||
double water_cut;
|
||||
if (std::abs(liquid_rate) != 0.) {
|
||||
water_cut = water_rate / liquid_rate;
|
||||
} else {
|
||||
water_cut = 0.0;
|
||||
}
|
||||
|
||||
const double max_water_cut_limit = econ_production_limits.maxWaterCut();
|
||||
if (water_cut > max_water_cut_limit) {
|
||||
water_cut_limit_violated = true;
|
||||
}
|
||||
|
||||
if (water_cut_limit_violated) {
|
||||
// need to handle the worst_offending_connection
|
||||
const int perf_start = map_entry[1];
|
||||
const int perf_number = map_entry[2];
|
||||
|
||||
std::vector<double> water_cut_perf(perf_number);
|
||||
for (int perf = 0; perf < perf_number; ++perf) {
|
||||
const int i_perf = perf_start + perf;
|
||||
const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
|
||||
const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
|
||||
const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
|
||||
if (std::abs(liquid_perf_rate) != 0.) {
|
||||
water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
|
||||
} else {
|
||||
water_cut_perf[perf] = 0.;
|
||||
}
|
||||
}
|
||||
|
||||
last_connection = (perf_number == 1);
|
||||
if (last_connection) {
|
||||
worst_offending_connection = 0;
|
||||
violation_extent = water_cut_perf[0] / max_water_cut_limit;
|
||||
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
|
||||
}
|
||||
|
||||
double max_water_cut_perf = 0.;
|
||||
for (int perf = 0; perf < perf_number; ++perf) {
|
||||
if (water_cut_perf[perf] > max_water_cut_perf) {
|
||||
worst_offending_connection = perf;
|
||||
max_water_cut_perf = water_cut_perf[perf];
|
||||
}
|
||||
}
|
||||
|
||||
assert(max_water_cut_perf != 0.);
|
||||
assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
|
||||
|
||||
violation_extent = max_water_cut_perf / max_water_cut_limit;
|
||||
}
|
||||
|
||||
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template <class WellState>
|
||||
void updateWellStateWithTarget(const WellControls* wc,
|
||||
const int current,
|
||||
const int well_index,
|
||||
WellState& xw) const
|
||||
{
|
||||
// number of phases
|
||||
const int np = wells().number_of_phases;
|
||||
// Updating well state and primary variables.
|
||||
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
|
||||
const double target = well_controls_iget_target(wc, current);
|
||||
const double* distr = well_controls_iget_distr(wc, current);
|
||||
switch (well_controls_iget_type(wc, current)) {
|
||||
case BHP:
|
||||
xw.bhp()[well_index] = target;
|
||||
break;
|
||||
|
||||
case THP: {
|
||||
double aqua = 0.0;
|
||||
double liquid = 0.0;
|
||||
double vapour = 0.0;
|
||||
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
|
||||
if (active_[ Water ]) {
|
||||
aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
vapour = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
|
||||
}
|
||||
|
||||
const int vfp = well_controls_iget_vfp(wc, current);
|
||||
const double& thp = well_controls_iget_target(wc, current);
|
||||
const double& alq = well_controls_iget_alq(wc, current);
|
||||
|
||||
//Set *BHP* target by calculating bhp from THP
|
||||
const WellType& well_type = wells().type[well_index];
|
||||
|
||||
// pick the density in the top layer
|
||||
const int perf = wells().well_connpos[well_index];
|
||||
const double rho = well_perforation_densities_[perf];
|
||||
|
||||
if (well_type == INJECTOR) {
|
||||
double dp = wellhelpers::computeHydrostaticCorrection(
|
||||
wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
|
||||
rho, gravity_);
|
||||
|
||||
xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
|
||||
}
|
||||
else if (well_type == PRODUCER) {
|
||||
double dp = wellhelpers::computeHydrostaticCorrection(
|
||||
wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
|
||||
rho, gravity_);
|
||||
|
||||
xw.bhp()[well_index] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
|
||||
}
|
||||
else {
|
||||
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
case RESERVOIR_RATE:
|
||||
// No direct change to any observable quantity at
|
||||
// surface condition. In this case, use existing
|
||||
// flow rates as initial conditions as reservoir
|
||||
// rate acts only in aggregate.
|
||||
break;
|
||||
|
||||
case SURFACE_RATE:
|
||||
// assign target value as initial guess for injectors and
|
||||
// single phase producers (orat, grat, wrat)
|
||||
const WellType& well_type = wells().type[well_index];
|
||||
if (well_type == INJECTOR) {
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
const double& compi = wells().comp_frac[np * well_index + phase];
|
||||
// TODO: it was commented out from the master branch already.
|
||||
//if (compi > 0.0) {
|
||||
xw.wellRates()[np*well_index + phase] = target * compi;
|
||||
//}
|
||||
}
|
||||
} else if (well_type == PRODUCER) {
|
||||
// only set target as initial rates for single phase
|
||||
// producers. (orat, grat and wrat, and not lrat)
|
||||
// lrat will result in numPhasesWithTargetsUnderThisControl == 2
|
||||
int numPhasesWithTargetsUnderThisControl = 0;
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
if (distr[phase] > 0.0) {
|
||||
numPhasesWithTargetsUnderThisControl += 1;
|
||||
}
|
||||
}
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
|
||||
xw.wellRates()[np*well_index + phase] = target * distr[phase];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
|
||||
}
|
||||
|
||||
break;
|
||||
} // end of switch
|
||||
|
||||
|
||||
std::vector<double> g = {1.0, 1.0, 0.01};
|
||||
if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
g[phase] = distr[phase];
|
||||
}
|
||||
}
|
||||
|
||||
// the number of wells
|
||||
const int nw = wells().number_of_wells;
|
||||
|
||||
switch (well_controls_iget_type(wc, current)) {
|
||||
case THP:
|
||||
case BHP: {
|
||||
const WellType& well_type = wells().type[well_index];
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] = 0.0;
|
||||
if (well_type == INJECTOR) {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] += xw.wellRates()[np*well_index + p] * wells().comp_frac[np*well_index + p];
|
||||
}
|
||||
} else {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] += g[p] * xw.wellRates()[np*well_index + p];
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
case RESERVOIR_RATE: // Intentional fall-through
|
||||
case SURFACE_RATE:
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] = xw.bhp()[well_index];
|
||||
break;
|
||||
} // end of switch
|
||||
|
||||
double tot_well_rate = 0.0;
|
||||
for (int p = 0; p < np; ++p) {
|
||||
tot_well_rate += g[p] * xw.wellRates()[np*well_index + p];
|
||||
}
|
||||
if(std::abs(tot_well_rate) > 0) {
|
||||
if (active_[ Water ]) {
|
||||
xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate;
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
|
||||
}
|
||||
} else {
|
||||
if (active_[ Water ]) {
|
||||
xw.wellSolutions()[WFrac*nw + well_index] = wells().comp_frac[np*well_index + Water];
|
||||
}
|
||||
|
||||
if (active_[ Gas ]) {
|
||||
xw.wellSolutions()[GFrac*nw + well_index] = wells().comp_frac[np*well_index + Gas];
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
WellState& xw) const;
|
||||
|
||||
};
|
||||
|
||||
|
||||
@ -504,7 +504,7 @@ namespace Opm {
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
resetWellControlFromState(WellState xw)
|
||||
resetWellControlFromState(WellState xw) const
|
||||
{
|
||||
const int nw = wells_->number_of_wells;
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
@ -617,6 +617,9 @@ namespace Opm {
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
typename StandardWellsDense<FluidSystem, BlackoilIndices>::EvalWell
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
@ -630,6 +633,9 @@ namespace Opm {
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
@ -646,6 +652,10 @@ namespace Opm {
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
@ -657,6 +667,10 @@ namespace Opm {
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
@ -1666,7 +1680,691 @@ namespace Opm {
|
||||
well_state.wellPotentials()[perf * np + p] = well_potentials[p].value();
|
||||
}
|
||||
}
|
||||
} // for (int w = 0; w < nw; ++w)
|
||||
} // end of for (int w = 0; w < nw; ++w)
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
WellCollection*
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
wellCollection() const
|
||||
{
|
||||
return well_collection_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
const std::vector<double>&
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
wellPerfEfficiencyFactors() const
|
||||
{
|
||||
return well_perforation_efficiency_factors_;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
calculateEfficiencyFactors()
|
||||
{
|
||||
if ( !localWellsActive() ) {
|
||||
return;
|
||||
}
|
||||
|
||||
const int nw = wells().number_of_wells;
|
||||
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
const std::string well_name = wells().name[w];
|
||||
const WellNode& well_node = wellCollection()->findWellNode(well_name);
|
||||
|
||||
const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor();
|
||||
|
||||
// assign the efficiency factor to each perforation related.
|
||||
for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w + 1]; ++perf) {
|
||||
well_perforation_efficiency_factors_[perf] = well_efficiency_factor;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
computeWellVoidageRates(const WellState& well_state,
|
||||
std::vector<double>& well_voidage_rates,
|
||||
std::vector<double>& voidage_conversion_coeffs) const
|
||||
{
|
||||
if ( !localWellsActive() ) {
|
||||
return;
|
||||
}
|
||||
// TODO: for now, we store the voidage rates for all the production wells.
|
||||
// For injection wells, the rates are stored as zero.
|
||||
// We only store the conversion coefficients for all the injection wells.
|
||||
// Later, more delicate model will be implemented here.
|
||||
// And for the moment, group control can only work for serial running.
|
||||
const int nw = well_state.numWells();
|
||||
const int np = well_state.numPhases();
|
||||
|
||||
// we calculate the voidage rate for each well, that means the sum of all the phases.
|
||||
well_voidage_rates.resize(nw, 0);
|
||||
// store the conversion coefficients, while only for the use of injection wells.
|
||||
voidage_conversion_coeffs.resize(nw * np, 1.0);
|
||||
|
||||
std::vector<double> well_rates(np, 0.0);
|
||||
std::vector<double> convert_coeff(np, 1.0);
|
||||
|
||||
for (int w = 0; w < nw; ++w) {
|
||||
const bool is_producer = wells().type[w] == PRODUCER;
|
||||
|
||||
// not sure necessary to change all the value to be positive
|
||||
if (is_producer) {
|
||||
std::transform(well_state.wellRates().begin() + np * w,
|
||||
well_state.wellRates().begin() + np * (w + 1),
|
||||
well_rates.begin(), std::negate<double>());
|
||||
|
||||
// the average hydrocarbon conditions of the whole field will be used
|
||||
const int fipreg = 0; // Not considering FIP for the moment.
|
||||
|
||||
rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff);
|
||||
well_voidage_rates[w] = std::inner_product(well_rates.begin(), well_rates.end(),
|
||||
convert_coeff.begin(), 0.0);
|
||||
} else {
|
||||
// TODO: Not sure whether will encounter situation with all zero rates
|
||||
// and whether it will cause problem here.
|
||||
std::copy(well_state.wellRates().begin() + np * w,
|
||||
well_state.wellRates().begin() + np * (w + 1),
|
||||
well_rates.begin());
|
||||
// the average hydrocarbon conditions of the whole field will be used
|
||||
const int fipreg = 0; // Not considering FIP for the moment.
|
||||
rate_converter_->calcCoeff(well_rates, fipreg, convert_coeff);
|
||||
std::copy(convert_coeff.begin(), convert_coeff.end(),
|
||||
voidage_conversion_coeffs.begin() + np * w);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
applyVREPGroupControl(WellState& well_state) const
|
||||
{
|
||||
if ( wellCollection()->havingVREPGroups() ) {
|
||||
std::vector<double> well_voidage_rates;
|
||||
std::vector<double> voidage_conversion_coeffs;
|
||||
computeWellVoidageRates(well_state, well_voidage_rates, voidage_conversion_coeffs);
|
||||
wellCollection()->applyVREPGroupControls(well_voidage_rates, voidage_conversion_coeffs);
|
||||
|
||||
// for the wells under group control, update the currentControls for the well_state
|
||||
for (const WellNode* well_node : wellCollection()->getLeafNodes()) {
|
||||
if (well_node->isInjector() && !well_node->individualControl()) {
|
||||
const int well_index = well_node->selfIndex();
|
||||
well_state.currentControls()[well_index] = well_node->groupControlIndex();
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
typename StandardWellsDense<FluidSystem, BlackoilIndices>::EvalWell
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
getBhp(const int wellIdx) const {
|
||||
const WellControls* wc = wells().ctrls[wellIdx];
|
||||
if (well_controls_get_current_type(wc) == BHP) {
|
||||
EvalWell bhp = 0.0;
|
||||
const double target_rate = well_controls_get_current_target(wc);
|
||||
bhp.setValue(target_rate);
|
||||
return bhp;
|
||||
} else if (well_controls_get_current_type(wc) == THP) {
|
||||
const int control = well_controls_get_current(wc);
|
||||
const double thp = well_controls_get_current_target(wc);
|
||||
const double alq = well_controls_iget_alq(wc, control);
|
||||
const int table_id = well_controls_iget_vfp(wc, control);
|
||||
EvalWell aqua = 0.0;
|
||||
EvalWell liquid = 0.0;
|
||||
EvalWell vapour = 0.0;
|
||||
EvalWell bhp = 0.0;
|
||||
double vfp_ref_depth = 0.0;
|
||||
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
|
||||
if (active_[ Water ]) {
|
||||
aqua = getQs(wellIdx, pu.phase_pos[ Water]);
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
liquid = getQs(wellIdx, pu.phase_pos[ Oil ]);
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
vapour = getQs(wellIdx, pu.phase_pos[ Gas ]);
|
||||
}
|
||||
if (wells().type[wellIdx] == INJECTOR) {
|
||||
bhp = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp);
|
||||
vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
|
||||
} else {
|
||||
bhp = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq);
|
||||
vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
|
||||
}
|
||||
|
||||
// pick the density in the top layer
|
||||
const int perf = wells().well_connpos[wellIdx];
|
||||
const double rho = well_perforation_densities_[perf];
|
||||
const double dp = wellhelpers::computeHydrostaticCorrection(wells(), wellIdx, vfp_ref_depth, rho, gravity_);
|
||||
bhp -= dp;
|
||||
return bhp;
|
||||
}
|
||||
|
||||
const int nw = wells().number_of_wells;
|
||||
return wellVariables_[nw*XvarWell + wellIdx];
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
typename StandardWellsDense<FluidSystem, BlackoilIndices>::EvalWell
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
getQs(const int wellIdx, const int phaseIdx) const
|
||||
{
|
||||
EvalWell qs = 0.0;
|
||||
const WellControls* wc = wells().ctrls[wellIdx];
|
||||
const int np = wells().number_of_phases;
|
||||
const int nw = wells().number_of_wells;
|
||||
const double target_rate = well_controls_get_current_target(wc);
|
||||
|
||||
if (wells().type[wellIdx] == INJECTOR) {
|
||||
const double comp_frac = wells().comp_frac[np*wellIdx + phaseIdx];
|
||||
if (comp_frac == 0.0) {
|
||||
return qs;
|
||||
}
|
||||
|
||||
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP) {
|
||||
return wellVariables_[nw*XvarWell + wellIdx];
|
||||
}
|
||||
qs.setValue(target_rate);
|
||||
return qs;
|
||||
}
|
||||
|
||||
// Producers
|
||||
if (well_controls_get_current_type(wc) == BHP || well_controls_get_current_type(wc) == THP ) {
|
||||
return wellVariables_[nw*XvarWell + wellIdx] * wellVolumeFractionScaled(wellIdx,phaseIdx);
|
||||
}
|
||||
|
||||
if (well_controls_get_current_type(wc) == SURFACE_RATE) {
|
||||
// checking how many phases are included in the rate control
|
||||
// to decide wheter it is a single phase rate control or not
|
||||
const double* distr = well_controls_get_current_distr(wc);
|
||||
int num_phases_under_rate_control = 0;
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
if (distr[phase] > 0.0) {
|
||||
num_phases_under_rate_control += 1;
|
||||
}
|
||||
}
|
||||
|
||||
// there should be at least one phase involved
|
||||
assert(num_phases_under_rate_control > 0);
|
||||
|
||||
// when it is a single phase rate limit
|
||||
if (num_phases_under_rate_control == 1) {
|
||||
if (distr[phaseIdx] == 1.0) {
|
||||
qs.setValue(target_rate);
|
||||
return qs;
|
||||
}
|
||||
|
||||
int currentControlIdx = 0;
|
||||
for (int i = 0; i < np; ++i) {
|
||||
currentControlIdx += wells().comp_frac[np*wellIdx + i] * i;
|
||||
}
|
||||
|
||||
const double eps = 1e-6;
|
||||
if (wellVolumeFractionScaled(wellIdx,currentControlIdx) < eps) {
|
||||
return qs;
|
||||
}
|
||||
return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / wellVolumeFractionScaled(wellIdx,currentControlIdx));
|
||||
}
|
||||
|
||||
// when it is a combined two phase rate limit, such like LRAT
|
||||
// we neec to calculate the rate for the certain phase
|
||||
if (num_phases_under_rate_control == 2) {
|
||||
EvalWell combined_volume_fraction = 0.;
|
||||
for (int p = 0; p < np; ++p) {
|
||||
if (distr[p] == 1.0) {
|
||||
combined_volume_fraction += wellVolumeFractionScaled(wellIdx, p);
|
||||
}
|
||||
}
|
||||
return (target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx) / combined_volume_fraction);
|
||||
}
|
||||
|
||||
// suppose three phase combined limit is the same with RESV
|
||||
// not tested yet.
|
||||
}
|
||||
|
||||
// ReservoirRate
|
||||
return target_rate * wellVolumeFractionScaled(wellIdx,phaseIdx);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
typename StandardWellsDense<FluidSystem, BlackoilIndices>::EvalWell
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
wellVolumeFraction(const int wellIdx, const int phaseIdx) const
|
||||
{
|
||||
const int nw = wells().number_of_wells;
|
||||
if (phaseIdx == Water) {
|
||||
return wellVariables_[WFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
if (phaseIdx == Gas) {
|
||||
return wellVariables_[GFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
// Oil fraction
|
||||
EvalWell well_fraction = 1.0;
|
||||
if (active_[Water]) {
|
||||
well_fraction -= wellVariables_[WFrac * nw + wellIdx];
|
||||
}
|
||||
|
||||
if (active_[Gas]) {
|
||||
well_fraction -= wellVariables_[GFrac * nw + wellIdx];
|
||||
}
|
||||
return well_fraction;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
typename StandardWellsDense<FluidSystem, BlackoilIndices>::EvalWell
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
wellVolumeFractionScaled(const int wellIdx, const int phaseIdx) const
|
||||
{
|
||||
const WellControls* wc = wells().ctrls[wellIdx];
|
||||
if (well_controls_get_current_type(wc) == RESERVOIR_RATE) {
|
||||
const double* distr = well_controls_get_current_distr(wc);
|
||||
return wellVolumeFraction(wellIdx, phaseIdx) / distr[phaseIdx];
|
||||
}
|
||||
std::vector<double> g = {1,1,0.01};
|
||||
return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
bool
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
checkRateEconLimits(const WellEconProductionLimits& econ_production_limits,
|
||||
const WellState& well_state,
|
||||
const int well_number) const
|
||||
{
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
const int np = well_state.numPhases();
|
||||
|
||||
if (econ_production_limits.onMinOilRate()) {
|
||||
assert(active_[Oil]);
|
||||
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
|
||||
const double min_oil_rate = econ_production_limits.minOilRate();
|
||||
if (std::abs(oil_rate) < min_oil_rate) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMinGasRate() ) {
|
||||
assert(active_[Gas]);
|
||||
const double gas_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Gas ] ];
|
||||
const double min_gas_rate = econ_production_limits.minGasRate();
|
||||
if (std::abs(gas_rate) < min_gas_rate) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMinLiquidRate() ) {
|
||||
assert(active_[Oil]);
|
||||
assert(active_[Water]);
|
||||
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
|
||||
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
|
||||
const double liquid_rate = oil_rate + water_rate;
|
||||
const double min_liquid_rate = econ_production_limits.minLiquidRate();
|
||||
if (std::abs(liquid_rate) < min_liquid_rate) {
|
||||
return true;
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMinReservoirFluidRate()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MIN_RESERVOIR_FLUID_RATE", "Minimum reservoir fluid production rate limit is not supported yet");
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
typename StandardWellsDense<FluidSystem, BlackoilIndices>::RatioCheckTuple
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
checkRatioEconLimits(const WellEconProductionLimits& econ_production_limits,
|
||||
const WellState& well_state,
|
||||
const WellMapEntryType& map_entry) const
|
||||
{
|
||||
// TODO: not sure how to define the worst-offending connection when more than one
|
||||
// ratio related limit is violated.
|
||||
// The defintion used here is that we define the violation extent based on the
|
||||
// ratio between the value and the corresponding limit.
|
||||
// For each violated limit, we decide the worst-offending connection separately.
|
||||
// Among the worst-offending connections, we use the one has the biggest violation
|
||||
// extent.
|
||||
|
||||
bool any_limit_violated = false;
|
||||
bool last_connection = false;
|
||||
int worst_offending_connection = INVALIDCONNECTION;
|
||||
double violation_extent = -1.0;
|
||||
|
||||
if (econ_production_limits.onMaxWaterCut()) {
|
||||
const RatioCheckTuple water_cut_return = checkMaxWaterCutLimit(econ_production_limits, well_state, map_entry);
|
||||
bool water_cut_violated = std::get<0>(water_cut_return);
|
||||
if (water_cut_violated) {
|
||||
any_limit_violated = true;
|
||||
const double violation_extent_water_cut = std::get<3>(water_cut_return);
|
||||
if (violation_extent_water_cut > violation_extent) {
|
||||
violation_extent = violation_extent_water_cut;
|
||||
worst_offending_connection = std::get<2>(water_cut_return);
|
||||
last_connection = std::get<1>(water_cut_return);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMaxGasOilRatio()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MAX_GOR", "the support for max Gas-Oil ratio is not implemented yet!");
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMaxWaterGasRatio()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MAX_WGR", "the support for max Water-Gas ratio is not implemented yet!");
|
||||
}
|
||||
|
||||
if (econ_production_limits.onMaxGasLiquidRatio()) {
|
||||
OpmLog::warning("NOT_SUPPORTING_MAX_GLR", "the support for max Gas-Liquid ratio is not implemented yet!");
|
||||
}
|
||||
|
||||
if (any_limit_violated) {
|
||||
assert(worst_offending_connection >=0);
|
||||
assert(violation_extent > 1.);
|
||||
}
|
||||
|
||||
return std::make_tuple(any_limit_violated, last_connection, worst_offending_connection, violation_extent);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
typename StandardWellsDense<FluidSystem, BlackoilIndices>::RatioCheckTuple
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
|
||||
const WellState& well_state,
|
||||
const WellMapEntryType& map_entry) const
|
||||
{
|
||||
bool water_cut_limit_violated = false;
|
||||
int worst_offending_connection = INVALIDCONNECTION;
|
||||
bool last_connection = false;
|
||||
double violation_extent = -1.0;
|
||||
|
||||
const int np = well_state.numPhases();
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
const int well_number = map_entry[0];
|
||||
|
||||
assert(active_[Oil]);
|
||||
assert(active_[Water]);
|
||||
|
||||
const double oil_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Oil ] ];
|
||||
const double water_rate = well_state.wellRates()[well_number * np + pu.phase_pos[ Water ] ];
|
||||
const double liquid_rate = oil_rate + water_rate;
|
||||
double water_cut;
|
||||
if (std::abs(liquid_rate) != 0.) {
|
||||
water_cut = water_rate / liquid_rate;
|
||||
} else {
|
||||
water_cut = 0.0;
|
||||
}
|
||||
|
||||
const double max_water_cut_limit = econ_production_limits.maxWaterCut();
|
||||
if (water_cut > max_water_cut_limit) {
|
||||
water_cut_limit_violated = true;
|
||||
}
|
||||
|
||||
if (water_cut_limit_violated) {
|
||||
// need to handle the worst_offending_connection
|
||||
const int perf_start = map_entry[1];
|
||||
const int perf_number = map_entry[2];
|
||||
|
||||
std::vector<double> water_cut_perf(perf_number);
|
||||
for (int perf = 0; perf < perf_number; ++perf) {
|
||||
const int i_perf = perf_start + perf;
|
||||
const double oil_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Oil ] ];
|
||||
const double water_perf_rate = well_state.perfPhaseRates()[i_perf * np + pu.phase_pos[ Water ] ];
|
||||
const double liquid_perf_rate = oil_perf_rate + water_perf_rate;
|
||||
if (std::abs(liquid_perf_rate) != 0.) {
|
||||
water_cut_perf[perf] = water_perf_rate / liquid_perf_rate;
|
||||
} else {
|
||||
water_cut_perf[perf] = 0.;
|
||||
}
|
||||
}
|
||||
|
||||
last_connection = (perf_number == 1);
|
||||
if (last_connection) {
|
||||
worst_offending_connection = 0;
|
||||
violation_extent = water_cut_perf[0] / max_water_cut_limit;
|
||||
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
|
||||
}
|
||||
|
||||
double max_water_cut_perf = 0.;
|
||||
for (int perf = 0; perf < perf_number; ++perf) {
|
||||
if (water_cut_perf[perf] > max_water_cut_perf) {
|
||||
worst_offending_connection = perf;
|
||||
max_water_cut_perf = water_cut_perf[perf];
|
||||
}
|
||||
}
|
||||
|
||||
assert(max_water_cut_perf != 0.);
|
||||
assert((worst_offending_connection >= 0) && (worst_offending_connection < perf_number));
|
||||
|
||||
violation_extent = max_water_cut_perf / max_water_cut_limit;
|
||||
}
|
||||
|
||||
return std::make_tuple(water_cut_limit_violated, last_connection, worst_offending_connection, violation_extent);
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
template<typename FluidSystem, typename BlackoilIndices>
|
||||
void
|
||||
StandardWellsDense<FluidSystem, BlackoilIndices>::
|
||||
updateWellStateWithTarget(const WellControls* wc,
|
||||
const int current,
|
||||
const int well_index,
|
||||
WellState& xw) const
|
||||
{
|
||||
// number of phases
|
||||
const int np = wells().number_of_phases;
|
||||
// Updating well state and primary variables.
|
||||
// Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
|
||||
const double target = well_controls_iget_target(wc, current);
|
||||
const double* distr = well_controls_iget_distr(wc, current);
|
||||
switch (well_controls_iget_type(wc, current)) {
|
||||
case BHP:
|
||||
xw.bhp()[well_index] = target;
|
||||
break;
|
||||
|
||||
case THP: {
|
||||
double aqua = 0.0;
|
||||
double liquid = 0.0;
|
||||
double vapour = 0.0;
|
||||
|
||||
const Opm::PhaseUsage& pu = phase_usage_;
|
||||
|
||||
if (active_[ Water ]) {
|
||||
aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
|
||||
}
|
||||
if (active_[ Oil ]) {
|
||||
liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
vapour = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
|
||||
}
|
||||
|
||||
const int vfp = well_controls_iget_vfp(wc, current);
|
||||
const double& thp = well_controls_iget_target(wc, current);
|
||||
const double& alq = well_controls_iget_alq(wc, current);
|
||||
|
||||
//Set *BHP* target by calculating bhp from THP
|
||||
const WellType& well_type = wells().type[well_index];
|
||||
|
||||
// pick the density in the top layer
|
||||
const int perf = wells().well_connpos[well_index];
|
||||
const double rho = well_perforation_densities_[perf];
|
||||
|
||||
if (well_type == INJECTOR) {
|
||||
const double dp = wellhelpers::computeHydrostaticCorrection(
|
||||
wells(), well_index, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
|
||||
rho, gravity_);
|
||||
|
||||
xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
|
||||
}
|
||||
else if (well_type == PRODUCER) {
|
||||
const double dp = wellhelpers::computeHydrostaticCorrection(
|
||||
wells(), well_index, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
|
||||
rho, gravity_);
|
||||
|
||||
xw.bhp()[well_index] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
|
||||
}
|
||||
else {
|
||||
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
case RESERVOIR_RATE:
|
||||
// No direct change to any observable quantity at
|
||||
// surface condition. In this case, use existing
|
||||
// flow rates as initial conditions as reservoir
|
||||
// rate acts only in aggregate.
|
||||
break;
|
||||
|
||||
case SURFACE_RATE:
|
||||
// assign target value as initial guess for injectors and
|
||||
// single phase producers (orat, grat, wrat)
|
||||
const WellType& well_type = wells().type[well_index];
|
||||
if (well_type == INJECTOR) {
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
const double& compi = wells().comp_frac[np * well_index + phase];
|
||||
// TODO: it was commented out from the master branch already.
|
||||
//if (compi > 0.0) {
|
||||
xw.wellRates()[np*well_index + phase] = target * compi;
|
||||
//}
|
||||
}
|
||||
} else if (well_type == PRODUCER) {
|
||||
// only set target as initial rates for single phase
|
||||
// producers. (orat, grat and wrat, and not lrat)
|
||||
// lrat will result in numPhasesWithTargetsUnderThisControl == 2
|
||||
int numPhasesWithTargetsUnderThisControl = 0;
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
if (distr[phase] > 0.0) {
|
||||
numPhasesWithTargetsUnderThisControl += 1;
|
||||
}
|
||||
}
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) {
|
||||
xw.wellRates()[np*well_index + phase] = target * distr[phase];
|
||||
}
|
||||
}
|
||||
} else {
|
||||
OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
|
||||
}
|
||||
|
||||
break;
|
||||
} // end of switch
|
||||
|
||||
|
||||
std::vector<double> g = {1.0, 1.0, 0.01};
|
||||
if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
|
||||
for (int phase = 0; phase < np; ++phase) {
|
||||
g[phase] = distr[phase];
|
||||
}
|
||||
}
|
||||
|
||||
// the number of wells
|
||||
const int nw = wells().number_of_wells;
|
||||
|
||||
switch (well_controls_iget_type(wc, current)) {
|
||||
case THP:
|
||||
case BHP: {
|
||||
const WellType& well_type = wells().type[well_index];
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] = 0.0;
|
||||
if (well_type == INJECTOR) {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] += xw.wellRates()[np*well_index + p] * wells().comp_frac[np*well_index + p];
|
||||
}
|
||||
} else {
|
||||
for (int p = 0; p < np; ++p) {
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] += g[p] * xw.wellRates()[np*well_index + p];
|
||||
}
|
||||
}
|
||||
break;
|
||||
}
|
||||
case RESERVOIR_RATE: // Intentional fall-through
|
||||
case SURFACE_RATE:
|
||||
xw.wellSolutions()[nw*XvarWell + well_index] = xw.bhp()[well_index];
|
||||
break;
|
||||
} // end of switch
|
||||
|
||||
double tot_well_rate = 0.0;
|
||||
for (int p = 0; p < np; ++p) {
|
||||
tot_well_rate += g[p] * xw.wellRates()[np*well_index + p];
|
||||
}
|
||||
if(std::abs(tot_well_rate) > 0) {
|
||||
if (active_[ Water ]) {
|
||||
xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate;
|
||||
}
|
||||
if (active_[ Gas ]) {
|
||||
xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
|
||||
}
|
||||
} else {
|
||||
if (active_[ Water ]) {
|
||||
xw.wellSolutions()[WFrac*nw + well_index] = wells().comp_frac[np*well_index + Water];
|
||||
}
|
||||
|
||||
if (active_[ Gas ]) {
|
||||
xw.wellSolutions()[GFrac*nw + well_index] = wells().comp_frac[np*well_index + Gas];
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
Loading…
Reference in New Issue
Block a user