diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp index 0b23077be..9cff12a49 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.cpp @@ -23,7 +23,6 @@ #include #include -#include #include #include @@ -38,10 +37,11 @@ #include #include -#include - #include +#include +#include + namespace Opm { template @@ -159,8 +159,7 @@ updateNewton(const BVectorWell& dwells, const double relaxation_factor, const double dFLimit, const bool stop_or_zero_rate_target, - const double max_pressure_change, - DeferredLogger& deferred_logger) + const double max_pressure_change) { const std::vector> old_primary_variables = value_; @@ -186,21 +185,8 @@ updateNewton(const BVectorWell& dwells, const double dx_limited = sign * std::min(std::abs(dwells[seg][SPres]) * relaxation_factor, max_pressure_change); // some cases might have defaulted bhp constraint of 1 bar, we use a slightly smaller value as the bhp lower limit for Newton update // so that bhp constaint can be an active control when needed. - if (seg == 0) { // top segment pressure is the bhp - constexpr double bhp_lower_limit = 1. * unit::barsa - 1. * unit::Pascal; - value_[seg][SPres] = std::max(old_primary_variables[seg][SPres] - dx_limited, bhp_lower_limit); - } else { - // enforcing 0. as the lower limit for the segment pressure of other segment - const double newton_update_value = old_primary_variables[seg][SPres] - dx_limited; - if (newton_update_value >= 0.) { - value_[seg][SPres] = newton_update_value; - } else { - const std::string msg = fmt::format("Segment with location index {} of well {} trying to acheive negative segment pressure {} bar" - " during the newton update, it will be chopped to 0.", seg, this->well_.name(), newton_update_value / unit::barsa); - deferred_logger.debug(msg); - value_[seg][SPres] = std::max(newton_update_value, 0.); - } - } + const double lower_limit = (seg == 0) ? bhp_lower_limit : seg_pres_lower_limit; + value_[seg][SPres] = std::max(old_primary_variables[seg][SPres] - dx_limited, lower_limit); } // update the total rate // TODO: should we have a limitation of the total rate change? @@ -642,6 +628,27 @@ getWQTotal() const return evaluation_[0][WQTotal]; } +template +void +MultisegmentWellPrimaryVariables:: +outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const +{ + std::string msg = fmt::format("outputting the segments for well {} with pressures close to the lower limits " + "for debugging purpose \n", this->well_.name()); + bool any_seg_pressure_close_to_limit = false; + for (std::size_t seg = 0; seg < value_.size(); ++seg) { + const double lower_limit = (seg == 0) ? bhp_lower_limit : seg_pres_lower_limit; + const double pres = Opm::getValue(this->getSegmentPressure(seg)); + if (pres <= std::numeric_limits::epsilon() + lower_limit) { + any_seg_pressure_close_to_limit = true; + fmt::format_to(std::back_inserter(msg), "seg {} : pressure {}\n", seg, pres / unit::barsa); + } + } + if (any_seg_pressure_close_to_limit) { + deferred_logger.debug(msg); + } +} + #define INSTANCE(...) \ template class MultisegmentWellPrimaryVariables,__VA_ARGS__,double>; diff --git a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp index e0077a91d..7fecb31a3 100644 --- a/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp +++ b/opm/simulators/wells/MultisegmentWellPrimaryVariables.hpp @@ -26,6 +26,7 @@ #include #include +#include #include #include @@ -97,8 +98,7 @@ public: const double relaxation_factor, const double DFLimit, const bool stop_or_zero_rate_target, - const double max_pressure_change, - DeferredLogger& deferred_logger); + const double max_pressure_change); //! \brief Copy values to well state. void copyToWellState(const MultisegmentWellGeneric& mswell, @@ -151,6 +151,9 @@ public: void setValue(const int idx, const std::array& val) { value_[idx] = val; } + //! output the segments with pressure close to lower pressure limit for debugging purpose + void outputLowLimitPressureSegments(DeferredLogger& deferred_logger) const; + private: //! \brief Handle non-reasonable fractions due to numerical overshoot. void processFractions(const int seg); @@ -168,6 +171,9 @@ private: std::vector> evaluation_; const WellInterfaceIndices& well_; //!< Reference to well interface + + static constexpr double bhp_lower_limit = 1. * unit::barsa - 1. * unit::Pascal; + static constexpr double seg_pres_lower_limit = 0.; }; } diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp index e257442cc..435de365b 100644 --- a/opm/simulators/wells/MultisegmentWell_impl.hpp +++ b/opm/simulators/wells/MultisegmentWell_impl.hpp @@ -689,8 +689,7 @@ namespace Opm relaxation_factor, dFLimit, stop_or_zero_rate_target, - max_pressure_change, - deferred_logger); + max_pressure_change); this->primary_variables_.copyToWellState(*this, getRefDensity(), stop_or_zero_rate_target, well_state, summary_state, deferred_logger); @@ -1762,6 +1761,7 @@ namespace Opm const std::string message = fmt::format(" Well {} did not converge in {} inner iterations (" "{} control/status switches).", this->name(), it, switch_count); deferred_logger.debug(message); + this->primary_variables_.outputLowLimitPressureSegments(deferred_logger); } return converged;