Report the operative CNV tolerance, the convergence gate and the oscillation source

opm/simulators/flow/ExtraConvergenceOutputThread.cpp

@@ -71,6 +71,23 @@ namespace {
         };
     }
 
+    // The trailing diagnostic columns, each with the longest token it can emit.
+    // Sizing them by their headings alone is not enough -- "GROUP+NETWORK" is
+    // wider than "Gate".  "WellStatus" appends free-form text after its token,
+    // which is why it must stay last.
+    auto diagnosticHeaders() noexcept
+    {
+        using namespace std::literals::string_view_literals;
+
+        return std::array {
+            std::pair { "CnvRelax"sv,   "FINALIT"sv  },  // to_string(CnvRelaxSource)
+            std::pair { "CnvTolUsed"sv, ""sv         },  // a number, never wider
+            std::pair { "OscSource"sv,  "WELLCTRL"sv },  // to_string(OscillationSource)
+            std::pair { "Gate"sv,       "GROUP+NETWORK"sv },
+            std::pair { "WellStatus"sv, "CONV"sv     },
+        };
+    }
+
     template <typename HeaderSequence>
     int maxHeaderSize(const HeaderSequence& headers)
     {
@@ -126,17 +143,27 @@ namespace {
             }
         }
 
-        const auto& metrics    = firstRequest.reports.front().reservoirConvergence();
-        const auto  headerSize = maxColHeaderSize(minColSize, getPhaseName, metrics);
+        const auto& metrics = firstRequest.reports.front().reservoirConvergence();
+
+        auto headerSize = maxColHeaderSize(minColSize, getPhaseName, metrics);
+        for (const auto& [heading, longestValue] : diagnosticHeaders()) {
+            headerSize = std::max({ headerSize,
+                                    static_cast<int>(heading.size() + 1),
+                                    static_cast<int>(longestValue.size() + 1) });
+        }
 
         for (const auto& metric : metrics) {
             os << std::right << std::setw(headerSize)
                << formatMetricColumn(getPhaseName, metric);
         }
 
+        for (const auto& [heading, _] : diagnosticHeaders()) {
+            os << std::right << std::setw(headerSize) << heading;
+        }
+
         // Note: Newline character intentionally placed in separate output
         // request to not influence right-justification of column header.
-        os << std::right << std::setw(headerSize) << "WellStatus" << '\n';
+        os << '\n';
 
         return { minColSize, headerSize };
     }
@@ -224,6 +251,46 @@ namespace {
         }
     }
 
+    void writeCnvRelaxation(std::ostream&                 os,
+                            const int                     colSize,
+                            const Opm::ConvergenceReport& report)
+    {
+        // Which condition granted the relaxed CNV tolerance, and the tolerance actually
+        // applied.  Without this the accepted state's quality is invisible: the same run
+        // may accept some steps at tolerance-cnv and others at its relaxed twin.
+        os << std::right << std::setw(colSize) << to_string(report.cnvRelaxSource());
+
+        const auto tol = report.cnvToleranceApplied();
+        if (tol > 0.0) {
+            os << std::right << std::setw(colSize) << tol;
+        }
+        else {
+            os << std::right << std::setw(colSize) << '-';
+        }
+    }
+
+    void writeOscillationSource(std::ostream&                 os,
+                                const int                     colSize,
+                                const Opm::ConvergenceReport& report)
+    {
+        os << std::right << std::setw(colSize) << to_string(report.oscillationSource());
+    }
+
+    void writeConvergenceGate(std::ostream&                 os,
+                              const int                     colSize,
+                              const Opm::ConvergenceReport& report)
+    {
+        // Which secondary gate, if any, kept this iteration from counting as
+        // converged.  These gates fire when the residual metrics themselves
+        // are already satisfied, so they identify iterations spent purely on
+        // well/group control changes or network balancing.
+        const auto gate = report.wellGroupTargetsViolated()
+            ? (report.networkNeedsMoreBalancing() ? "GROUP+NETWORK" : "GROUP")
+            : (report.networkNeedsMoreBalancing() ? "NETWORK" : "NONE");
+
+        os << std::right << std::setw(colSize) << gate;
+    }
+
     void writeConvergenceRequest(std::ostream&                                           os,
                                  const Opm::ConvergenceOutputThread::ConvertToTimeUnits& convertTime,
                                  const int                                               firstColSize,
@@ -243,6 +310,9 @@ namespace {
             writePenaltyCount(os, firstColSize, report);
 
             writeReservoirConvergence(os, colSize, report);
+            writeCnvRelaxation(os, colSize, report);
+            writeOscillationSource(os, colSize, report);
+            writeConvergenceGate(os, colSize, report);
             writeWellConvergence(os, colSize, report);
 
             os << '\n';

opm/simulators/flow/NonlinearSystemBlackOilReservoir_impl.hpp

@@ -173,6 +173,11 @@ initialLinearization(SimulatorReportSingle& report,
         auto convrep = getConvergence(timer, maxIter, residual_norms);
         report.converged = convrep.converged() &&
                            this->simulator_.problem().iterationContext().iteration() >= minIter;
+        if (report.converged &&
+            convrep.cnvRelaxSource() != ConvergenceReport::CnvRelaxSource::None)
+        {
+            ++report.relaxed_cnv_acceptances;
+        }
         ConvergenceReport::Severity severity = convrep.severityOfWorstFailure();
         this->convergence_reports_.back().report.push_back(std::move(convrep));
 
@@ -290,8 +295,21 @@ nonlinearIterationNewton(const SimulatorTimerInterface& timer,
             if (isOscillate) {
                 this->current_relaxation_ -= nonlinear_solver.relaxIncrement();
                 this->current_relaxation_ = std::max(this->current_relaxation_, nonlinear_solver.relaxMax());
+                // The detector reads reservoir residual history and the response damps the
+                // reservoir update, but the oscillation may originate in the well/group
+                // control layer.  Record what was unsatisfied here so the two can be told
+                // apart.
+                auto source = ConvergenceReport::OscillationSource::NotDetected;
+                if (!this->convergence_reports_.empty() &&
+                    !this->convergence_reports_.back().report.empty())
+                {
+                    auto& convrep = this->convergence_reports_.back().report.back();
+                    source = classifyOscillationSource(convrep);
+                    convrep.setOscillationSource(source);
+                }
                 if (this->terminalOutputEnabled()) {
-                    OpmLog::info("    Oscillating behavior detected: Relaxation set to "
+                    OpmLog::info("    Oscillating behavior detected (" + to_string(source)
+                                 + "): Relaxation set to "
                                  + std::to_string(this->current_relaxation_));
                 }
             }
@@ -816,6 +834,18 @@ getReservoirConvergence(const double reportTime,
 
     const auto tol_cnv = use_relaxed_cnv ? tolerance_cnv_relaxed : this->param_.tolerance_cnv_;
     const auto tol_mb  = use_relaxed_mb ? this->param_.tolerance_mb_relaxed_ : this->param_.tolerance_mb_;
+
+    // Record which condition granted the relaxation, so the accepted state's quality is
+    // visible in the output. Order mirrors the code's own precedence.
+    {
+        using RS = ConvergenceReport::CnvRelaxSource;
+        const auto source = relax_dsol_cnv          ? RS::SolChange
+                          : relax_pv_fraction_cnv   ? RS::PvFraction
+                          : relax_final_iteration_cnv ? RS::FinalIter
+                          : relax_iter_cnv          ? RS::IterCount
+                                                    : RS::None;
+        report.setCnvRelaxation(source, static_cast<double>(tol_cnv));
+    }
     const auto tol_cnv_energy = use_relaxed_cnv ? this->param_.tolerance_cnv_energy_relaxed_ : this->param_.tolerance_cnv_energy_;
     const auto tol_eb = use_relaxed_mb ?  this->param_.tolerance_energy_balance_relaxed_ : this->param_.tolerance_energy_balance_;
 

opm/simulators/timestepping/ConvergenceReport.cpp

@@ -19,6 +19,7 @@
 
 #include <opm/simulators/timestepping/ConvergenceReport.hpp>
 
+#include <algorithm>
 #include <cassert>
 #include <stdexcept>
 #include <string>
@@ -107,4 +108,57 @@ namespace Opm
                            pc.nonConverged, pc.distanceDecay, pc.largeWellResiduals, pc.total());
     }
 
+    std::string to_string(const ConvergenceReport::CnvRelaxSource s)
+    {
+        using S = ConvergenceReport::CnvRelaxSource;
+        switch (s) {
+        case S::None:       return "NONE";
+        case S::PvFraction: return "PVFRAC";
+        case S::SolChange:  return "DSOL";
+        case S::FinalIter:  return "FINALIT";
+        case S::IterCount:  return "ITER";
+        }
+        return "?";
+    }
+
+    std::string to_string(const ConvergenceReport::OscillationSource s)
+    {
+        using S = ConvergenceReport::OscillationSource;
+        switch (s) {
+        case S::NotDetected: return "NONE";
+        case S::Reservoir:   return "RES";
+        case S::WellControl: return "WELLCTRL";
+        case S::Mixed:       return "MIXED";
+        }
+        return "?";
+    }
+
+    ConvergenceReport::OscillationSource
+    classifyOscillationSource(const ConvergenceReport& report)
+    {
+        using WT = ConvergenceReport::WellFailure::Type;
+
+        const auto& wellFailures = report.wellFailures();
+        const bool controlUnsatisfied =
+            std::any_of(wellFailures.begin(), wellFailures.end(),
+                        [](const ConvergenceReport::WellFailure& wf)
+                        {
+                            return (wf.type() == WT::ControlBHP)
+                                || (wf.type() == WT::ControlTHP)
+                                || (wf.type() == WT::ControlRate);
+                        })
+            || report.wellGroupTargetsViolated()
+            || report.networkNeedsMoreBalancing();
+
+        const bool reservoirUnsatisfied = report.reservoirFailed();
+
+        if (controlUnsatisfied && reservoirUnsatisfied) {
+            return ConvergenceReport::OscillationSource::Mixed;
+        }
+        if (controlUnsatisfied) {
+            return ConvergenceReport::OscillationSource::WellControl;
+        }
+        return ConvergenceReport::OscillationSource::Reservoir;
+    }
+
 } // namespace Opm

opm/simulators/timestepping/ConvergenceReport.hpp

@@ -240,6 +240,29 @@ namespace Opm
             std::string well_name_ {};
         };
 
+        /// Which condition, if any, granted the relaxed CNV tolerance this iteration.
+        /// Ordered by precedence in the code so a single value can be reported even
+        /// though several conditions may hold at once.
+        enum struct CnvRelaxSource {
+            None,        //!< strict tolerance applied
+            PvFraction,  //!< non-strict pore volume below relaxed_max_pv_fraction
+            SolChange,   //!< TUNINGDP-style solution-change targets met (disables CNV)
+            FinalIter,   //!< min_strict_cnv_iter < 0 and iteration == NewtonMaxIterations
+            IterCount,   //!< min_strict_cnv_iter >= 0 and iteration >= it
+        };
+
+        /// What was unsatisfied at the iteration where the Newton solver detected an
+        /// oscillation.  An association, not a proof of cause.  Note the base rate when
+        /// reading this: a non-final iteration almost always has unsatisfied reservoir
+        /// residuals, so Reservoir is the common label and the discriminating signal is
+        /// how often the control layer *also* participates (Mixed + WellControl).
+        enum struct OscillationSource {
+            NotDetected,
+            Reservoir,    //!< reservoir residuals only
+            WellControl,  //!< well control equations or a group/network gate only
+            Mixed,        //!< both
+        };
+
         // ----------- Mutating member functions -----------
 
         ConvergenceReport()
@@ -305,6 +328,21 @@ namespace Opm
             this->eligiblePoreVolume_ = eligiblePoreVolume;
         }
 
+        //! \brief Record that the Newton solver damped this iteration for oscillation,
+        //! together with what was unsatisfied at that point.
+        void setOscillationSource(const OscillationSource source)
+        {
+            this->oscillationSource_ = source;
+        }
+
+        //! \brief Record which condition granted the relaxed CNV tolerance and the
+        //! tolerance actually applied, so an accepted step's quality is visible.
+        void setCnvRelaxation(const CnvRelaxSource source, const double toleranceApplied)
+        {
+            this->cnvRelaxSource_ = source;
+            this->cnvToleranceApplied_ = toleranceApplied;
+        }
+
         ConvergenceReport& operator+=(const ConvergenceReport& other)
         {
             reportTime_ = std::max(reportTime_, other.reportTime_);
@@ -331,6 +369,13 @@ namespace Opm
                 this->eligiblePoreVolume_ = other.eligiblePoreVolume_;
             }
 
+            // Same reasoning as the pore-volume split: only the reservoir report sets
+            // these, so take 'other's values whenever it carries them.
+            if (other.cnvToleranceApplied_ > 0.0) {
+                this->cnvRelaxSource_ = other.cnvRelaxSource_;
+                this->cnvToleranceApplied_ = other.cnvToleranceApplied_;
+            }
+
             return *this;
         }
 
@@ -351,13 +396,46 @@ namespace Opm
             return this->cnvPvSplit_;
         }
 
+        //! \brief What was unsatisfied when oscillation damping was triggered, if it was.
+        OscillationSource oscillationSource() const
+        {
+            return this->oscillationSource_;
+        }
+
+        //! \brief Which condition granted the relaxed CNV tolerance this iteration.
+        CnvRelaxSource cnvRelaxSource() const
+        {
+            return this->cnvRelaxSource_;
+        }
+
+        //! \brief The CNV tolerance actually applied (0 if not recorded).
+        double cnvToleranceApplied() const
+        {
+            return this->cnvToleranceApplied_;
+        }
+
         bool converged() const
         {
             return (status_ == AllGood)
                 && !wellGroupTargetsViolated_
                 && !network_needs_more_balancing_force_another_newton_iteration_;
         }
 
+        //! \brief Whether well/group control targets were violated (a control
+        //! or target changed this iteration), forcing another Newton iteration
+        //! even when all residual metrics are converged.
+        bool wellGroupTargetsViolated() const
+        {
+            return wellGroupTargetsViolated_;
+        }
+
+        //! \brief Whether the network balance forced another Newton iteration
+        //! even when all residual metrics are converged.
+        bool networkNeedsMoreBalancing() const
+        {
+            return network_needs_more_balancing_force_another_newton_iteration_;
+        }
+
         bool reservoirFailed() const
         {
             return status_ & ReservoirFailed;
@@ -437,6 +515,9 @@ namespace Opm
             serializer(this->network_needs_more_balancing_force_another_newton_iteration_);
             serializer(this->cnvPvSplit_);
             serializer(this->eligiblePoreVolume_);
+            serializer(this->cnvRelaxSource_);
+            serializer(this->cnvToleranceApplied_);
+            serializer(this->oscillationSource_);
             serializer(this->penaltyCard_);
         }
 
@@ -454,6 +535,9 @@ namespace Opm
         bool network_needs_more_balancing_force_another_newton_iteration_;
         CnvPvSplit cnvPvSplit_{};
         double eligiblePoreVolume_{};
+        CnvRelaxSource cnvRelaxSource_{CnvRelaxSource::None};
+        double cnvToleranceApplied_{};
+        OscillationSource oscillationSource_{OscillationSource::NotDetected};
         PenaltyCard penaltyCard_;
     };
 
@@ -474,6 +558,13 @@ namespace Opm
 
     std::string to_string(const ConvergenceReport::PenaltyCard& pc);
 
+    std::string to_string(const ConvergenceReport::CnvRelaxSource s);
+
+    std::string to_string(const ConvergenceReport::OscillationSource s);
+
+    /// Classify what was unsatisfied in a report, for oscillation-source reporting.
+    ConvergenceReport::OscillationSource classifyOscillationSource(const ConvergenceReport& report);
+
 
 
 } // namespace Opm