diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index f4ebac8f6..d181e186a 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -608,77 +608,71 @@ namespace Opm checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits, const WellState& well_state) const { - bool water_cut_limit_violated = false; - int worst_offending_completion = INVALIDCOMPLETION; - double violation_extent = -1.0; - - const int np = number_of_phases_; - const Opm::PhaseUsage& pu = phaseUsage(); - const int well_number = index_of_well_; assert(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)); assert(FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)); - 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; - } + auto waterCut = [](const double oil_rate, const double water_rate) { + + // both rate should be in the same direction + assert(oil_rate * water_rate >= 0.); + + const double liquid_rate = oil_rate + water_rate; + if (liquid_rate != 0.) { + return (water_rate / liquid_rate); + } else { + return 0.; + } + }; + + const int np = number_of_phases_; + const Opm::PhaseUsage& pu = phaseUsage(); + + const double oil_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Oil ] ]; + const double water_rate = well_state.wellRates()[index_of_well_ * np + pu.phase_pos[ Water ] ]; + const double water_cut = waterCut(oil_rate, water_rate); const double max_water_cut_limit = econ_production_limits.maxWaterCut(); - if (water_cut > max_water_cut_limit) { - water_cut_limit_violated = true; - } + assert(max_water_cut_limit != 0.); + + const bool water_cut_limit_violated = (water_cut > max_water_cut_limit); + + int worst_offending_completion = INVALIDCOMPLETION; + double violation_extent = -1.0; if (water_cut_limit_violated) { - // need to handle the worst_offending_connection - const int perf_start = first_perf_; - const int perf_number = number_of_perforations_; + // the maximum water cut value of the completions + double max_water_cut_completion = 0.; - 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.; + // look for the worst_offending_completion + for (const auto& completion : completions_) { + + double oil_completion_rate = 0.; + double water_completion_rate = 0.; + + const std::vector& conns = completion.second; + for(const int c : conns) { + const int index_con = c + first_perf_; + + const double oil_connection_rate = well_state.perfPhaseRates()[index_con * np + pu.phase_pos[ Oil ] ]; + oil_completion_rate += oil_connection_rate; + + const double water_connection_rate = well_state.perfPhaseRates()[index_con * np + pu.phase_pos[ Water ] ]; + water_completion_rate += water_connection_rate; } - } - const auto& completions = well_ecl_.getCompletions(); - const auto& connections = well_ecl_.getConnections(); - int complnumIdx = 0; - std::vector water_cut_in_completions(completions.size(), 0.0); - for (const auto& completion : completions) { - int complnum = completion.first; - for (int perf = 0; perf < perf_number; ++perf) { - if (complnum == connections.get ( perf ).complnum()) { - water_cut_in_completions[complnumIdx] += water_cut_perf[perf]; - } - } - complnumIdx++; - } + const double water_cut_completion = waterCut(oil_completion_rate, water_completion_rate); - double max_water_cut_perf = 0.; - complnumIdx = 0; - for (const auto& completion : completions) { - if (water_cut_in_completions[complnumIdx] > max_water_cut_perf) { + if (water_cut_completion > max_water_cut_completion) { worst_offending_completion = completion.first; - max_water_cut_perf = water_cut_in_completions[complnumIdx]; + max_water_cut_completion = water_cut_completion; } - complnumIdx++; - } + } // end of for (const auto& completion : completions_) - assert(max_water_cut_limit != 0.); + assert(max_water_cut_completion >= max_water_cut_limit); assert(worst_offending_completion != INVALIDCOMPLETION); - violation_extent = max_water_cut_perf / max_water_cut_limit; + + violation_extent = max_water_cut_completion / max_water_cut_limit; } return std::make_tuple(water_cut_limit_violated, worst_offending_completion, violation_extent);