From 5f344eef266eb2da2fbb7b9ae58cc4727e1cf95e Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 20 Jun 2023 23:14:50 +0200 Subject: [PATCH] avoid dividing by zero in updateWellStateRates() --- opm/simulators/wells/WellInterface_impl.hpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 0d0c6482f..5167795de 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1181,11 +1181,13 @@ namespace Opm // Set the currently-zero phase flows to be nonzero in proportion to well_q_s. const double initial_nonzero_rate = ws.surface_rates[nonzero_rate_index]; const int comp_idx_nz = this->flowPhaseToEbosCompIdx(nonzero_rate_index); - for (int p = 0; p < this->number_of_phases_; ++p) { - if (p != nonzero_rate_index) { - const int comp_idx = this->flowPhaseToEbosCompIdx(p); - double& rate = ws.surface_rates[p]; - rate = (initial_nonzero_rate/well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]); + if (std::abs(well_q_s[comp_idx_nz]) > floating_point_error_epsilon) { + for (int p = 0; p < this->number_of_phases_; ++p) { + if (p != nonzero_rate_index) { + const int comp_idx = this->flowPhaseToEbosCompIdx(p); + double& rate = ws.surface_rates[p]; + rate = (initial_nonzero_rate / well_q_s[comp_idx_nz]) * (well_q_s[comp_idx]); + } } } }