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 @@ class BlackoilModelNldd {
                     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

@@ -94,6 +94,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>
@@ -228,6 +232,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

@@ -137,6 +137,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 {
@@ -284,6 +289,13 @@ struct BlackoilModelParameters
 
     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