diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp index b49c9428f..8a7368977 100644 --- a/opm/autodiff/StandardWellsDense.hpp +++ b/opm/autodiff/StandardWellsDense.hpp @@ -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& - 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& well_voidage_rates, - std::vector& 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 well_rates(np, 0.0); - std::vector 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()); - - // 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 well_voidage_rates; - std::vector 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& 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 g = {1,1,0.01}; - return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]); - } - - - - template 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 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 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 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 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 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; }; diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp index b78bf8236..8c6aae79d 100644 --- a/opm/autodiff/StandardWellsDense_impl.hpp +++ b/opm/autodiff/StandardWellsDense_impl.hpp @@ -504,7 +504,7 @@ namespace Opm { template void StandardWellsDense:: - 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 StandardWellsDense::EvalWell StandardWellsDense:: @@ -630,6 +633,9 @@ namespace Opm { } + + + template void StandardWellsDense:: @@ -646,6 +652,10 @@ namespace Opm { } } + + + + template void StandardWellsDense:: @@ -657,6 +667,10 @@ namespace Opm { } } + + + + template void StandardWellsDense:: @@ -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 + WellCollection* + StandardWellsDense:: + wellCollection() const + { + return well_collection_; + } + + + + + template + const std::vector& + StandardWellsDense:: + wellPerfEfficiencyFactors() const + { + return well_perforation_efficiency_factors_; + } + + + + + + template + void + StandardWellsDense:: + 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 + void + StandardWellsDense:: + computeWellVoidageRates(const WellState& well_state, + std::vector& well_voidage_rates, + std::vector& 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 well_rates(np, 0.0); + std::vector 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()); + + // 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 + void + StandardWellsDense:: + applyVREPGroupControl(WellState& well_state) const + { + if ( wellCollection()->havingVREPGroups() ) { + std::vector well_voidage_rates; + std::vector 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 StandardWellsDense::EvalWell + StandardWellsDense:: + 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 StandardWellsDense::EvalWell + StandardWellsDense:: + 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 StandardWellsDense::EvalWell + StandardWellsDense:: + 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 StandardWellsDense::EvalWell + StandardWellsDense:: + 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 g = {1,1,0.01}; + return (wellVolumeFraction(wellIdx, phaseIdx) / g[phaseIdx]); + } + + + + + + template + bool + StandardWellsDense:: + 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 StandardWellsDense::RatioCheckTuple + StandardWellsDense:: + 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 StandardWellsDense::RatioCheckTuple + StandardWellsDense:: + 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 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 + void + StandardWellsDense:: + 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 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