diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.cpp b/opm/simulators/wells/StandardWellPrimaryVariables.cpp index ef0385adc..4d387ecb0 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.cpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.cpp @@ -40,19 +40,36 @@ #include #include +#include + #include #include +#include +#include namespace { - constexpr double EPSILON = 1.0e-14; - -//! \brief Relaxation factor considering only one fraction value. +/** + * @brief Relaxation factor considering only one fraction value. + * @param old_value The previous value before the change. + * @param dx The solved newton update for the value. + * @param well_name The name of the well that calls this function + * @param value_name The name/meaning of the old_value. + * @param deferred_logger The deferred logger used for logging results and messages. + * @return The calculated relaxation factor. + */ template Scalar relaxationFactorFraction(const Scalar old_value, - const Scalar dx) + const Scalar dx, + const std::string& well_name, + const std::string& value_name, + Opm::DeferredLogger& deferred_logger) { - assert(old_value >= -EPSILON && old_value <= (1.0+EPSILON)); + constexpr double epislon = std::numeric_limits::epsilon(); + if (old_value < -epislon || old_value > 1.0 + epislon) { + const std::string msg = fmt::format(" illegal fraction value {} {} is found for well {}", value_name, old_value, well_name); + OPM_DEFLOG_THROW(Opm::NumericalProblem, msg, deferred_logger); + } const Scalar& safe_old_value = std::clamp(old_value, 0.0, 1.0); Scalar relaxation_factor = 1.; @@ -228,15 +245,17 @@ template void StandardWellPrimaryVariables:: updateNewton(const BVectorWell& dwells, const bool stop_or_zero_rate_target, - [[maybe_unused]] const double dFLimit, - const double dBHPLimit) + const double dFLimit, + const double dBHPLimit, + DeferredLogger& deferred_logger) { // for injectors, very typical one of the fractions will be one, and it is easy to get zero value // fractions. not sure what is the best way to handle it yet, so we just use 1.0 here // The relaxationFactorFractionProducer code does not take into account solvent // so we use 1.0 for cases with solvent. [[maybe_unused]] const double relaxation_factor_fractions = - (well_.isProducer() && !Indices::enableSolvent) ? this->relaxationFactorFractionsProducer(dwells) : 1.0; + (well_.isProducer() && !Indices::enableSolvent) ? this->relaxationFactorFractionsProducer(dwells, deferred_logger) + : 1.0; // update the second and third well variable (The flux fractions) if constexpr (has_wfrac_variable) { @@ -642,7 +661,7 @@ processFractions() template double StandardWellPrimaryVariables:: -relaxationFactorFractionsProducer(const BVectorWell& dwells) const +relaxationFactorFractionsProducer(const BVectorWell& dwells, DeferredLogger& deferred_logger) const { // TODO: not considering solvent yet // 0.95 is a experimental value, which remains to be optimized @@ -651,13 +670,19 @@ relaxationFactorFractionsProducer(const BVectorWell& dwells) const if (FluidSystem::numActivePhases() > 1) { if constexpr (has_wfrac_variable) { const double relaxation_factor_w = relaxationFactorFraction(value_[WFrac], - dwells[0][WFrac]); + dwells[0][WFrac], + this->well_.name(), + "WFrac", + deferred_logger); relaxation_factor = std::min(relaxation_factor, relaxation_factor_w); } if constexpr (has_gfrac_variable) { const double relaxation_factor_g = relaxationFactorFraction(value_[GFrac], - dwells[0][GFrac]); + dwells[0][GFrac], + this->well_.name(), + "GFrac", + deferred_logger); relaxation_factor = std::min(relaxation_factor, relaxation_factor_g); } @@ -670,7 +695,7 @@ relaxationFactorFractionsProducer(const BVectorWell& dwells) const const double possible_updated_sum = original_sum - relaxed_update; // We only relax if fraction is above 1. // The newton solver should handle the rest - const double epsilon = 0.001; + constexpr double epsilon = 0.001; if (possible_updated_sum > 1.0 + epsilon) { // since the orignal sum <= 1.0 the epsilon asserts that // the relaxed_update is non trivial. @@ -680,7 +705,10 @@ relaxationFactorFractionsProducer(const BVectorWell& dwells) const relaxation_factor *= further_relaxation_factor; } } - assert(relaxation_factor >= 0.0 && relaxation_factor <= 1.0); + if (relaxation_factor < 0.0 || relaxation_factor > 1.0) { + const std::string msg = fmt::format(" illegal relaxation factor {} is obtained for well {}", relaxation_factor, this->well_.name()); + OPM_DEFLOG_THROW(NumericalProblem, msg, deferred_logger); + } } return relaxation_factor; } diff --git a/opm/simulators/wells/StandardWellPrimaryVariables.hpp b/opm/simulators/wells/StandardWellPrimaryVariables.hpp index 10037f677..4f4f7de09 100644 --- a/opm/simulators/wells/StandardWellPrimaryVariables.hpp +++ b/opm/simulators/wells/StandardWellPrimaryVariables.hpp @@ -112,7 +112,8 @@ public: void updateNewton(const BVectorWell& dwells, const bool stop_or_zero_rate_target, const double dFLimit, - const double dBHPLimit); + const double dBHPLimit, + DeferredLogger& deferred_logger); //! \brief Update polymer molecular weight values from newton update vector. void updateNewtonPolyMW(const BVectorWell& dwells); @@ -150,7 +151,7 @@ public: private: //! \brief Calculate a relaxation factor for producers. //! \details To avoid overshoot of the fractions which might result in negative rates. - double relaxationFactorFractionsProducer(const BVectorWell& dwells) const; + double relaxationFactorFractionsProducer(const BVectorWell& dwells, DeferredLogger& deferred_logger) const; //! \brief Returns volume fraction for a component. EvalWell volumeFraction(const unsigned compIdx) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index a34f54b53..d7c753a7d 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -703,7 +703,7 @@ namespace Opm { const double dFLimit = this->param_.dwell_fraction_max_; const double dBHPLimit = this->param_.dbhp_max_rel_; - this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit); + this->primary_variables_.updateNewton(dwells, stop_or_zero_rate_target, dFLimit, dBHPLimit, deferred_logger); // for the water velocity and skin pressure if constexpr (Base::has_polymermw) {