From a828ce336e1513810dda04d41980651d8f1af09f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 29 Apr 2025 09:55:50 +0200 Subject: [PATCH] Ensure that the initialization is not producing above-bhp-limit rates. --- opm/simulators/wells/WellInterface_impl.hpp | 27 ++++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/opm/simulators/wells/WellInterface_impl.hpp b/opm/simulators/wells/WellInterface_impl.hpp index 13b0ba0dfb..7aefc36342 100644 --- a/opm/simulators/wells/WellInterface_impl.hpp +++ b/opm/simulators/wells/WellInterface_impl.hpp @@ -1650,25 +1650,40 @@ namespace Opm } if (nonzero_rate_index == -1) { - // No nonzero rates. + // No nonzero rates on input. // Use the computed rate directly, or scaled by a factor // 0.5 (to avoid too high values) if it was evaluated at 1 bar. - const double factor = rates_evaluated_at_1bar ? 0.5 : 1.0; + const Scalar factor = rates_evaluated_at_1bar ? 0.5 : 1.0; for (int p = 0; p < this->number_of_phases_; ++p) { ws.surface_rates[p] = factor * well_q_s[p]; } return; } - // Set the currently-zero phase flows to be nonzero in proportion to well_q_s. + // If we are here, we had a single nonzero rate for the well, + // typically from a rate constraint. We must make sure it is + // respected, so if it was lower than the calculated rate for + // the same phase we scale all rates to match. const Scalar initial_nonzero_rate = ws.surface_rates[nonzero_rate_index]; - if (std::abs(well_q_s[nonzero_rate_index]) > floating_point_error_epsilon) { + const Scalar computed_rate = well_q_s[nonzero_rate_index]; + if (std::abs(initial_nonzero_rate) < std::abs(computed_rate)) { + // Note that both rates below are negative. The factor should be < 1.0. + const Scalar factor = initial_nonzero_rate / computed_rate; + assert(factor < 1.0); for (int p = 0; p < this->number_of_phases_; ++p) { + // We skip the nonzero_rate_index, as that should remain as it was. if (p != nonzero_rate_index) { - Scalar& rate = ws.surface_rates[p]; - rate = (initial_nonzero_rate / well_q_s[nonzero_rate_index]) * (well_q_s[p]); + ws.surface_rates[p] = factor * well_q_s[p]; } } + return; + } + + // If we are here, we had a single nonzero rate, but it was + // higher than the one calculated from the bhp limit, so we + // use the calculated rates. + for (int p = 0; p < this->number_of_phases_; ++p) { + ws.surface_rates[p] = well_q_s[p]; } }