Merge pull request #5590 from jakobtorben/convergence_monitors

Convergence monitors

opm/simulators/flow/BlackoilModel.hpp

@@ -366,9 +366,14 @@ namespace Opm {
 
                 // Throw if any NaN or too large residual found.
                 if (severity == ConvergenceReport::Severity::NotANumber) {
+                    failureReport_ += report;
                     OPM_THROW_PROBLEM(NumericalProblem, "NaN residual found!");
                 } else if (severity == ConvergenceReport::Severity::TooLarge) {
+                    failureReport_ += report;
                     OPM_THROW_NOLOG(NumericalProblem, "Too large residual found!");
+                } else if (severity == ConvergenceReport::Severity::ConvergenceMonitorFailure) {
+                    failureReport_ += report;
+                    OPM_THROW_PROBLEM(ConvergenceMonitorFailure, "Total penalty count exceeded cut-off-limit of " + std::to_string(param_.convergence_monitoring_cutoff_));
                 }
             }
             report.update_time += perfTimer.stop();
@@ -392,6 +397,9 @@ namespace Opm {
                 // For each iteration we store in a vector the norms of the residual of
                 // the mass balance for each active phase, the well flux and the well equations.
                 residual_norms_history_.clear();
+                total_penaltyCard_.reset();
+                prev_above_tolerance_ = 0;
+                prev_distance_ = std::numeric_limits<double>::infinity();
                 current_relaxation_ = 1.0;
                 dx_old_ = 0.0;
                 convergence_reports_.push_back({timer.reportStepNum(), timer.currentStepNum(), {}});
@@ -1062,7 +1070,7 @@ namespace Opm {
                         report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
                     }
 
-                    report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
+                    report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
                 }
             }
 
@@ -1108,6 +1116,51 @@ namespace Opm {
             return report;
         }
 
+    void checkCardPenalty(ConvergenceReport& report, int iteration)
+    {
+
+        const auto& current_metrics = report.reservoirConvergence();
+        auto distances = std::vector<double>(current_metrics.size(), 0.0);
+        int current_above_tolerance = 0;
+
+        for (size_t i = 0; i < current_metrics.size(); ++i) {
+            distances[i] = std::max(std::log10(current_metrics[i].value()/current_metrics[i].tolerance()), 0.0);
+                // Count number of metrics above tolerance
+                if (current_metrics[i].value() > current_metrics[i].tolerance()) {
+                    current_above_tolerance++;
+                }
+            }
+
+        // use L1 norm of the distances vector
+        double current_distance = std::accumulate(distances.begin(), distances.end(), 0.0);
+
+        if (iteration > 0) {
+            // Add penalty if number of metrics above tolerance has increased
+            if (current_above_tolerance > prev_above_tolerance_) {
+                report.addNonConvergedPenalty();
+            }
+
+            if (current_distance > param_.convergence_monitoring_decay_factor_ * prev_distance_) {
+                report.addDistanceDecayPenalty();
+            }
+        }
+
+        prev_distance_ = current_distance;
+        prev_above_tolerance_ = current_above_tolerance;
+
+        if (report.wellFailures().size() > 0) {
+            report.addLargeWellResidualsPenalty();
+        }
+
+        total_penaltyCard_ += report.getPenaltyCard();
+
+        if (param_.convergence_monitoring_ && (total_penaltyCard_.total() > param_.convergence_monitoring_cutoff_)) {
+            report.setReservoirFailed({ConvergenceReport::ReservoirFailure::Type::ConvergenceMonitorFailure,
+                                       ConvergenceReport::Severity::ConvergenceMonitorFailure,
+                                       -1}); // -1 indicates it's not specific to any component
+        }
+    }
+
         /// Compute convergence based on total mass balance (tol_mb) and maximum
         /// residual mass balance (tol_cnv).
         /// \param[in]   timer       simulation timer
@@ -1129,6 +1182,9 @@ namespace Opm {
                 OPM_TIMEBLOCK(getWellConvergence);
                 report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged());
             }
+
+            checkCardPenalty(report, iteration);
+
             return report;
         }
 
@@ -1395,7 +1451,9 @@ namespace Opm {
         Scalar drMaxRel() const { return param_.dr_max_rel_; }
         Scalar maxResidualAllowed() const { return param_.max_residual_allowed_; }
         double linear_solve_setup_time_;
-
+        ConvergenceReport::PenaltyCard total_penaltyCard_;
+        double prev_distance_ = std::numeric_limits<double>::infinity();
+        int prev_above_tolerance_ = 0;
     public:
         std::vector<bool> wasSwitched_;
     };

opm/simulators/flow/BlackoilModelNldd.hpp

@@ -678,7 +678,7 @@ private:
                     report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
                 }
 
-                report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
+                report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii], tol[ii]);
             }
         }
 

opm/simulators/flow/BlackoilModelParameters.cpp

@@ -95,6 +95,10 @@ BlackoilModelParameters<Scalar>::BlackoilModelParameters()
     network_max_iterations_ = Parameters::Get<Parameters::NetworkMaxIterations>();
     local_domain_ordering_ = domainOrderingMeasureFromString(Parameters::Get<Parameters::LocalDomainsOrderingMeasure>());
     write_partitions_ = Parameters::Get<Parameters::DebugEmitCellPartition>();
+
+    convergence_monitoring_ = Parameters::Get<Parameters::ConvergenceMonitoring>();
+    convergence_monitoring_cutoff_ = Parameters::Get<Parameters::ConvergenceMonitoringCutOff>();
+    convergence_monitoring_decay_factor_ = Parameters::Get<Parameters::ConvergenceMonitoringDecayFactor<Scalar>>();
 }
 
 template<class Scalar>
@@ -231,6 +235,13 @@ void BlackoilModelParameters<Scalar>::registerParameters()
     Parameters::Register<Parameters::DebugEmitCellPartition>
         ("Whether or not to emit cell partitions as a debugging aid.");
 
+    Parameters::Register<Parameters::ConvergenceMonitoring>
+        ("Enable convergence monitoring");
+    Parameters::Register<Parameters::ConvergenceMonitoringCutOff>
+        ("Cut off limit for convergence monitoring");
+    Parameters::Register<Parameters::ConvergenceMonitoringDecayFactor<Scalar>>
+        ("Decay factor for convergence monitoring");
+
     Parameters::Hide<Parameters::DebugEmitCellPartition>();
 
     // if openMP is available, determine the number threads per process automatically.

opm/simulators/flow/BlackoilModelParameters.hpp

@@ -138,6 +138,11 @@ struct LocalDomainsPartitioningImbalance { static constexpr Scalar value = 1.03;
 struct LocalDomainsPartitioningMethod { static constexpr auto value = "zoltan"; };
 struct LocalDomainsOrderingMeasure { static constexpr auto value = "maxpressure"; };
 
+struct ConvergenceMonitoring { static constexpr bool value = false; };
+struct ConvergenceMonitoringCutOff { static constexpr int value = 6; };
+template<class Scalar>
+struct ConvergenceMonitoringDecayFactor { static constexpr Scalar value = 0.75; };
+
 } // namespace Opm::Parameters
 
 namespace Opm {
@@ -288,6 +293,13 @@ public:
 
     bool write_partitions_{false};
 
+    /// Whether to enable convergence monitoring
+    bool convergence_monitoring_;
+    /// Cut-off limit for convergence monitoring
+    int convergence_monitoring_cutoff_;
+    /// Decay factor used in convergence monitoring
+    Scalar convergence_monitoring_decay_factor_;
+
     /// Construct from user parameters or defaults.
     BlackoilModelParameters();
 

opm/simulators/flow/ExtraConvergenceOutputThread.cpp

@@ -65,6 +65,9 @@ namespace {
             "CnvCellCntConv"s,
             "CnvCellCntRelax"s,
             "CnvCellCntUnconv"s,
+            "PenaltyNonConv"s,
+            "PenaltyDecay"s,
+            "PenaltyWellRes"s,
         };
     }
 
@@ -190,6 +193,17 @@ namespace {
         }
     }
 
+    void writePenaltyCount(std::ostream&                 os,
+                           const std::string::size_type  firstColSize,
+                           const Opm::ConvergenceReport& report)
+    {
+        const auto& penaltyCard = report.getPenaltyCard();
+
+        os << ' ' << std::setw(firstColSize) << penaltyCard.nonConverged;
+        os << ' ' << std::setw(firstColSize) << penaltyCard.distanceDecay;
+        os << ' ' << std::setw(firstColSize) << penaltyCard.largeWellResiduals;
+    }
+
     void writeReservoirConvergence(std::ostream&                 os,
                                    const std::string::size_type  colSize,
                                    const Opm::ConvergenceReport& report)
@@ -229,6 +243,8 @@ namespace {
 
             writeCnvPvSplit(os, expectNumCnvSplit, firstColSize, report);
 
+            writePenaltyCount(os, firstColSize, report);
+
             writeReservoirConvergence(os, colSize, report);
             writeWellConvergence(os, colSize, report);
 

opm/simulators/timestepping/AdaptiveTimeStepping.hpp

@@ -247,6 +247,10 @@ void registerAdaptiveParameters();
                     logException_(e, solverVerbose_);
                     // since linearIterations is < 0 this will restart the solver
                 }
+                catch (const ConvergenceMonitorFailure& e) {
+                    substepReport = solver.failureReport();
+                    causeOfFailure = "Convergence monitor failure";
+                }
                 catch (const LinearSolverProblem& e) {
                     substepReport = solver.failureReport();
                     causeOfFailure = "Linear solver convergence failure";

opm/simulators/timestepping/ConvergenceReport.cpp

@@ -59,6 +59,8 @@ namespace Opm
             return "TooLarge";
         case S::NotANumber:
             return "NotANumber";
+        case S::ConvergenceMonitorFailure:
+            return "ConvergenceMonitorFailure";
         }
         throw std::logic_error("Unknown ConvergenceReport::Severity");
     }
@@ -97,4 +99,10 @@ namespace Opm
     }
 
 
+    std::string to_string(const ConvergenceReport::PenaltyCard& pc)
+    {
+        return fmt::format("PenaltyCard {{ NonConverged: {}, DistanceDecay: {}, LargeWellResiduals: {}, Total: {} }}",
+                           pc.nonConverged, pc.distanceDecay, pc.largeWellResiduals, pc.total());
+    }
+
 } // namespace Opm

opm/simulators/timestepping/ConvergenceReport.hpp

@@ -45,12 +45,45 @@ namespace Opm
             ReservoirFailed = 1 << 0,
             WellFailed      = 1 << 1,
         };
-
+        // More severe problems should have higher numbers
         enum struct Severity {
             None       = 0,
             Normal     = 1,
-            TooLarge   = 2,
-            NotANumber = 3,
+            ConvergenceMonitorFailure = 2,
+            TooLarge   = 3,
+            NotANumber = 4,
+        };
+
+        struct PenaltyCard {
+            int nonConverged{0};
+            int distanceDecay{0};
+            int largeWellResiduals{0};
+
+            int total() const {
+                return nonConverged + distanceDecay + largeWellResiduals;
+            }
+
+            void reset()
+            {
+                nonConverged = 0;
+                distanceDecay = 0;
+                largeWellResiduals = 0;
+            }
+
+            PenaltyCard& operator+=(const PenaltyCard& other) {
+                nonConverged += other.nonConverged;
+                distanceDecay += other.distanceDecay;
+                largeWellResiduals += other.largeWellResiduals;
+                return *this;
+            }
+
+            template <typename Serializer>
+            void serializeOp(Serializer& serializer)
+            {
+                serializer(nonConverged);
+                serializer(distanceDecay);
+                serializer(largeWellResiduals);
+            }
         };
 
         using CnvPvSplit = std::pair<
@@ -60,7 +93,7 @@ namespace Opm
         class ReservoirFailure
         {
         public:
-            enum struct Type { Invalid, MassBalance, Cnv };
+            enum struct Type { Invalid, MassBalance, Cnv, ConvergenceMonitorFailure };
 
             // Default constructor needed for object serialisation.  Don't
             // use this for anything else.
@@ -97,13 +130,14 @@ namespace Opm
             // use this for anything else.
             ReservoirConvergenceMetric() = default;
 
-            ReservoirConvergenceMetric(ReservoirFailure::Type t, int phase, double value)
-                : type_(t), phase_(phase), value_(value)
+            ReservoirConvergenceMetric(ReservoirFailure::Type t, int phase, double value, double tolerance)
+                : type_(t), phase_(phase), value_(value), tolerance_(tolerance)
             {}
 
             ReservoirFailure::Type type() const { return type_; }
             int phase() const { return phase_; }
             double value() const { return value_; }
+            double tolerance() const { return tolerance_; }
 
             template <typename Serializer>
             void serializeOp(Serializer& serializer)
@@ -111,6 +145,7 @@ namespace Opm
                 serializer(this->type_);
                 serializer(this->phase_);
                 serializer(this->value_);
+                serializer(this->tolerance_);
             }
 
         private:
@@ -119,6 +154,7 @@ namespace Opm
             ReservoirFailure::Type type_ { ReservoirFailure::Type::Invalid };
             int phase_ { -1 };
             double value_ { 0.0 };
+            double tolerance_ { 0.0 };
         };
 
         class WellFailure
@@ -166,6 +202,43 @@ namespace Opm
             std::string well_name_ {};
         };
 
+        class WellConvergenceMetric
+        {
+        public:
+            // Default constructor needed for object serialisation.  Don't
+            // use this for anything else.
+            WellConvergenceMetric() = default;
+
+            WellConvergenceMetric(WellFailure::Type t, Severity s, int phase, double value, const std::string& well_name)
+                : type_(t), severity_(s), phase_(phase), value_(value), well_name_(well_name)
+            {}
+
+            WellFailure::Type type() const { return type_; }
+            Severity severity() const { return severity_; }
+            int phase() const { return phase_; }
+            double value() const { return value_; }
+            const std::string& wellName() const { return well_name_; }
+
+            template <typename Serializer>
+            void serializeOp(Serializer& serializer)
+            {
+                serializer(this->type_);
+                serializer(this->severity_);
+                serializer(this->phase_);
+                serializer(this->value_);
+                serializer(this->well_name_);
+            }
+
+        private:
+            // Note to maintainers: If you change this list of data members,
+            // then please update serializeOp() accordingly.
+            WellFailure::Type type_ { WellFailure::Type::Invalid };
+            Severity severity_ { Severity::None };
+            int phase_ { -1 };
+            double value_ { 0.0 };
+            std::string well_name_ {};
+        };
+
         // ----------- Mutating member functions -----------
 
         ConvergenceReport()
@@ -206,6 +279,12 @@ namespace Opm
             this->res_convergence_.emplace_back(std::forward<Args>(args)...);
         }
 
+        template <typename... Args>
+        void setWellConvergenceMetric(Args&&... args)
+        {
+            this->well_convergence_.emplace_back(std::forward<Args>(args)...);
+        }
+
         void setWellGroupTargetsViolated(const bool wellGroupTargetsViolated)
         {
             wellGroupTargetsViolated_ = wellGroupTargetsViolated;
@@ -225,6 +304,7 @@ namespace Opm
             res_failures_.insert(res_failures_.end(), other.res_failures_.begin(), other.res_failures_.end());
             well_failures_.insert(well_failures_.end(), other.well_failures_.begin(), other.well_failures_.end());
             res_convergence_.insert(res_convergence_.end(), other.res_convergence_.begin(), other.res_convergence_.end());
+            well_convergence_.insert(well_convergence_.end(), other.well_convergence_.begin(), other.well_convergence_.end());
             assert(reservoirFailed() != res_failures_.empty());
             assert(wellFailed() != well_failures_.empty());
             wellGroupTargetsViolated_ = (wellGroupTargetsViolated_ || other.wellGroupTargetsViolated_);
@@ -291,6 +371,31 @@ namespace Opm
             return well_failures_;
         }
 
+        const std::vector<WellConvergenceMetric>& wellConvergence() const
+        {
+            return well_convergence_;
+        }
+
+        const PenaltyCard& getPenaltyCard() const
+        {
+            return penaltyCard_;
+        }
+
+        void addNonConvergedPenalty()
+        {
+            penaltyCard_.nonConverged++;
+        }
+
+        void addDistanceDecayPenalty()
+        {
+            penaltyCard_.distanceDecay++;
+        }
+
+        void addLargeWellResidualsPenalty()
+        {
+            penaltyCard_.largeWellResiduals++;
+        }
+
         Severity severityOfWorstFailure() const
         {
             // A function to get the worst of two severities.
@@ -315,9 +420,11 @@ namespace Opm
             serializer(this->res_failures_);
             serializer(this->well_failures_);
             serializer(this->res_convergence_);
+            serializer(this->well_convergence_);
             serializer(this->wellGroupTargetsViolated_);
             serializer(this->cnvPvSplit_);
             serializer(this->eligiblePoreVolume_);
+            serializer(this->penaltyCard_);
         }
 
     private:
@@ -329,9 +436,11 @@ namespace Opm
         std::vector<ReservoirFailure> res_failures_;
         std::vector<WellFailure> well_failures_;
         std::vector<ReservoirConvergenceMetric> res_convergence_;
+        std::vector<WellConvergenceMetric> well_convergence_;
         bool wellGroupTargetsViolated_;
         CnvPvSplit cnvPvSplit_{};
         double eligiblePoreVolume_{};
+        PenaltyCard penaltyCard_;
     };
 
     struct StepReport
@@ -349,6 +458,9 @@ namespace Opm
 
     std::string to_string(const ConvergenceReport::WellFailure& wf);
 
+    std::string to_string(const ConvergenceReport::PenaltyCard& pc);
+
+
 
 } // namespace Opm
 

opm/simulators/wells/StandardWellEval.cpp

@@ -134,13 +134,20 @@ getWellConvergence(const WellState<Scalar>& well_state,
 
         if (std::isnan(well_flux_residual[compIdx])) {
             report.setWellFailed({type, CR::Severity::NotANumber, compIdx, baseif_.name()});
+            report.setWellConvergenceMetric(type, CR::Severity::NotANumber, compIdx, well_flux_residual[compIdx], baseif_.name());
         } else if (well_flux_residual[compIdx] > maxResidualAllowed) {
             report.setWellFailed({type, CR::Severity::TooLarge, compIdx, baseif_.name()});
+            report.setWellConvergenceMetric(type, CR::Severity::TooLarge, compIdx, well_flux_residual[compIdx], baseif_.name());
         } else if (!relax_tolerance && well_flux_residual[compIdx] > tol_wells) {
             report.setWellFailed({type, CR::Severity::Normal, compIdx, baseif_.name()});
+            report.setWellConvergenceMetric(type, CR::Severity::Normal, compIdx, well_flux_residual[compIdx], baseif_.name());
         } else if (well_flux_residual[compIdx] > relaxed_tolerance_flow) {
             report.setWellFailed({type, CR::Severity::Normal, compIdx, baseif_.name()});
+            report.setWellConvergenceMetric(type, CR::Severity::Normal, compIdx, well_flux_residual[compIdx], baseif_.name());
+        } else {
+            report.setWellConvergenceMetric(CR::WellFailure::Type::Invalid, CR::Severity::None, compIdx, well_flux_residual[compIdx], baseif_.name());
         }
+
     }
 
     WellConvergence(baseif_).