From edf8786a586df7fdfe15b16b0569b96b6d001fa9 Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Fri, 8 Apr 2022 11:07:34 +0200
Subject: [PATCH] avoid throw for nan residuals instead return well as
 unconverged

---
 opm/simulators/wells/MultisegmentWellEval.cpp | 16 +++++++-------
 opm/simulators/wells/MultisegmentWellEval.hpp |  5 +++--
 .../wells/MultisegmentWell_impl.hpp           | 21 ++++++++++++++++---
 3 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/opm/simulators/wells/MultisegmentWellEval.cpp b/opm/simulators/wells/MultisegmentWellEval.cpp
index 72469cfbb..419782fdf 100644
--- a/opm/simulators/wells/MultisegmentWellEval.cpp
+++ b/opm/simulators/wells/MultisegmentWellEval.cpp
@@ -1619,10 +1619,10 @@ assemblePressureEq(const int seg,
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
-std::vector<Scalar>
+std::pair<bool, std::vector<Scalar> >
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
-getWellResiduals(const std::vector<Scalar>& B_avg,
-                 DeferredLogger& deferred_logger) const
+getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
+                       DeferredLogger& deferred_logger) const
 {
     assert(int(B_avg.size() ) == baseif_.numComponents());
     std::vector<Scalar> residuals(numWellEq + 1, 0.0);
@@ -1638,8 +1638,9 @@ getWellResiduals(const std::vector<Scalar>& B_avg,
                 }
             }
             if (std::isnan(residual) || std::isinf(residual)) {
-                OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for residal get for well " << baseif_.name()
-                                                << " segment " << seg << " eq_idx " << eq_idx, deferred_logger);
+                deferred_logger.debug("nan or inf value for residal get for well " + baseif_.name()
+                                      + " segment " + std::to_string(seg) + " eq_idx " + std::to_string(eq_idx));
+                return {false, residuals};
             }
 
             if (residual > residuals[eq_idx]) {
@@ -1652,12 +1653,13 @@ getWellResiduals(const std::vector<Scalar>& B_avg,
     {
         const double control_residual = std::abs(resWell_[0][numWellEq - 1]);
         if (std::isnan(control_residual) || std::isinf(control_residual)) {
-           OPM_DEFLOG_THROW(NumericalIssue, "nan or inf value for control residal get for well " << baseif_.name(), deferred_logger);
+           deferred_logger.debug("nan or inf value for control residal get for well " + baseif_.name());
+           return {false, residuals};
         }
         residuals[numWellEq] = control_residual;
     }
 
-    return residuals;
+    return {true, residuals};
 }
 
 template<typename FluidSystem, typename Indices, typename Scalar>
diff --git a/opm/simulators/wells/MultisegmentWellEval.hpp b/opm/simulators/wells/MultisegmentWellEval.hpp
index 0654e1e77..aba4bf324 100644
--- a/opm/simulators/wells/MultisegmentWellEval.hpp
+++ b/opm/simulators/wells/MultisegmentWellEval.hpp
@@ -199,8 +199,9 @@ protected:
     EvalWell getWQTotal() const;
 
 
-    std::vector<Scalar> getWellResiduals(const std::vector<Scalar>& B_avg,
-                                         DeferredLogger& deferred_logger) const;
+    std::pair<bool, std::vector<Scalar> >
+    getFiniteWellResiduals(const std::vector<Scalar>& B_avg,
+                           DeferredLogger& deferred_logger) const;
 
     double getControlTolerance(const WellState& well_state,
                                const double tolerance_wells,
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index 112f44e24..51c8f5c05 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -1369,7 +1369,14 @@ namespace Opm
 
         const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
         const WellState well_state0 = well_state;
-        const std::vector<Scalar> residuals0 = this->getWellResiduals(Base::B_avg_, deferred_logger);
+
+        {
+            // getWellFiniteResiduals returns false for nan/inf residuals
+            const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
+            if(!isFinite)
+                return false;
+        }
+
         std::vector<std::vector<Scalar> > residual_history;
         std::vector<double> measure_history;
         int it = 0;
@@ -1394,12 +1401,20 @@ namespace Opm
                 break;
             }
 
-            residual_history.push_back(this->getWellResiduals(Base::B_avg_, deferred_logger));
-            measure_history.push_back(this->getResidualMeasureValue(well_state,
+            {
+                // getFinteWellResiduals returns false for nan/inf residuals
+                const auto& [isFinite, residuals] = this->getFiniteWellResiduals(Base::B_avg_, deferred_logger);
+                if (!isFinite)
+                    return false;
+
+                residual_history.push_back(residuals);
+                measure_history.push_back(this->getResidualMeasureValue(well_state,
                                                                     residual_history[it],
                                                                     this->param_.tolerance_wells_,
                                                                     this->param_.tolerance_pressure_ms_wells_,
                                                                     deferred_logger) );
+            }
+
 
             bool is_oscillate = false;
             bool is_stagnate = false;