Capture Timestep's Non-Linear Convergence History

This enables outputting a formatted record of the limiting MB and
CNV quantities as time and non-linear iterations progress.  This, in
turn, is intended for diagnostic and analysis purposes and will not
be output unless specifically requested.

In particular, add a new type,

    ConvergenceReport::ReservoirConvergenceMetric

which captures the convergence metric type (MB or CNV) along with
the associate phase and numerical value of the convergence metric.
We add a vector of these convergence metric objects as a new data
member of the ConvergenceReport.

Finally, foreshadowing the intended use case, also store the report
time in the ConvergenceReport object.

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -888,7 +888,8 @@ namespace Opm {
             return grid_.comm().sum(errorPV);
         }
 
-        ConvergenceReport getReservoirConvergence(const double dt,
+        ConvergenceReport getReservoirConvergence(const double reportTime,
+                                                  const double dt,
                                                   const int iteration,
                                                   std::vector<Scalar>& B_avg,
                                                   std::vector<Scalar>& residual_norms)
@@ -927,7 +928,7 @@ namespace Opm {
             }
 
             // Create convergence report.
-            ConvergenceReport report;
+            ConvergenceReport report{reportTime};
             using CR = ConvergenceReport;
             for (int compIdx = 0; compIdx < numComp; ++compIdx) {
                 double res[2] = { mass_balance_residual[compIdx], CNV[compIdx] };
@@ -953,6 +954,7 @@ namespace Opm {
                     } else if (res[ii] > tol[ii]) {
                         report.setReservoirFailed({types[ii], CR::Severity::Normal, compIdx});
                     }
+                    report.setReservoirConvergenceMetric(types[ii], compIdx, res[ii]);
                 }
             }
 
@@ -1003,7 +1005,9 @@ namespace Opm {
         {
             // Get convergence reports for reservoir and wells.
             std::vector<Scalar> B_avg(numEq, 0.0);
-            auto report = getReservoirConvergence(timer.currentStepLength(), iteration, B_avg, residual_norms);
+            auto report = getReservoirConvergence(timer.simulationTimeElapsed(),
+                                                  timer.currentStepLength(),
+                                                  iteration, B_avg, residual_norms);
             report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged());
 
             return report;

opm/simulators/timestepping/ConvergenceReport.hpp

@@ -21,9 +21,11 @@
 #ifndef OPM_CONVERGENCEREPORT_HEADER_INCLUDED
 #define OPM_CONVERGENCEREPORT_HEADER_INCLUDED
 
+#include <algorithm>
 #include <cassert>
 #include <numeric>
 #include <string>
+#include <utility>
 #include <vector>
 
 namespace Opm
@@ -61,6 +63,21 @@ namespace Opm
             Severity severity_;
             int phase_;
         };
+        class ReservoirConvergenceMetric
+        {
+        public:
+            ReservoirConvergenceMetric(ReservoirFailure::Type t, int phase, double value)
+                : type_(t), phase_(phase), value_(value)
+            {
+            }
+            ReservoirFailure::Type type() const { return type_; }
+            int phase() const { return phase_; }
+            double value() const { return value_; }
+        private:
+            ReservoirFailure::Type type_;
+            int phase_;
+            double value_;
+        };
         class WellFailure
         {
         public:
@@ -83,7 +100,13 @@ namespace Opm
         // ----------- Mutating member functions -----------
 
         ConvergenceReport()
-            : status_{AllGood}
+            : ConvergenceReport{0.0}
+        {
+        }
+
+        explicit ConvergenceReport(const double reportTime)
+            : reportTime_{reportTime}
+            , status_{AllGood}
             , res_failures_{}
             , well_failures_{}
             , wellGroupTargetsViolated_(false)
@@ -110,6 +133,12 @@ namespace Opm
             well_failures_.push_back(wf);
         }
 
+        template <typename... Args>
+        void setReservoirConvergenceMetric(Args&&... args)
+        {
+            this->res_convergence_.emplace_back(std::forward<Args>(args)...);
+        }
+
         void setWellGroupTargetsViolated(const bool wellGroupTargetsViolated)
         {
             wellGroupTargetsViolated_ = wellGroupTargetsViolated;
@@ -117,9 +146,11 @@ namespace Opm
 
         ConvergenceReport& operator+=(const ConvergenceReport& other)
         {
+            reportTime_ = std::max(reportTime_, other.reportTime_);
             status_ = static_cast<Status>(status_ | other.status_);
             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());
             assert(reservoirFailed() != res_failures_.empty());
             assert(wellFailed() != well_failures_.empty());
             wellGroupTargetsViolated_ = (wellGroupTargetsViolated_ || other.wellGroupTargetsViolated_);
@@ -128,6 +159,11 @@ namespace Opm
 
         // ----------- Const member functions (queries) -----------
 
+        double reportTime() const
+        {
+            return reportTime_;
+        }
+
         bool converged() const
         {
             return (status_ == AllGood) && !wellGroupTargetsViolated_;
@@ -148,6 +184,11 @@ namespace Opm
             return res_failures_;
         }
 
+        const std::vector<ReservoirConvergenceMetric>& reservoirConvergence() const
+        {
+            return res_convergence_;
+        }
+
         const std::vector<WellFailure>& wellFailures() const
         {
             return well_failures_;
@@ -172,9 +213,11 @@ namespace Opm
     private:
 
         // ----------- Member variables -----------
+        double reportTime_;
         Status status_;
         std::vector<ReservoirFailure> res_failures_;
         std::vector<WellFailure> well_failures_;
+        std::vector<ReservoirConvergenceMetric> res_convergence_;
         bool wellGroupTargetsViolated_;
     };
 

opm/simulators/timestepping/gatherConvergenceReport.cpp

@@ -1,6 +1,6 @@
 /*
+  Copyright 2018, 2022 Equinor ASA.
   Copyright 2018 SINTEF Digital, Mathematics and Cybernetics.
-  Copyright 2018 Equinor.
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -43,6 +43,18 @@ namespace
         MPI_Pack(&phase, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator);
     }
 
+    void packReservoirConvergenceMetric(const ConvergenceReport::ReservoirConvergenceMetric& m,
+                                        std::vector<char>& buf,
+                                        int& offset, MPI_Comm mpi_communicator)
+    {
+        int type = static_cast<int>(m.type());
+        int phase = m.phase();
+        double value = m.value();
+        MPI_Pack(&type, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator);
+        MPI_Pack(&phase, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator);
+        MPI_Pack(&value, 1, MPI_DOUBLE, buf.data(), buf.size(), &offset, mpi_communicator);
+    }
+
     void packWellFailure(const ConvergenceReport::WellFailure& f,
                          std::vector<char>& buf,
                          int& offset, MPI_Comm mpi_communicator)
@@ -65,12 +77,21 @@ namespace
         // Pack the data.
         // Status will not be packed, it is possible to deduce from the other data.
         // Reservoir failures.
+        double reportTime = local_report.reportTime();
+        MPI_Pack(&reportTime, 1, MPI_DOUBLE, buf.data(), buf.size(), &offset, mpi_communicator);
         const auto rf = local_report.reservoirFailures();
         int num_rf = rf.size();
         MPI_Pack(&num_rf, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator);
         for (const auto& f : rf) {
             packReservoirFailure(f, buf, offset, mpi_communicator);
         }
+        // Reservoir convergence metrics.
+        const auto rm = local_report.reservoirConvergence();
+        int num_rm = rm.size();
+        MPI_Pack(&num_rm, 1, MPI_INT, buf.data(), buf.size(), &offset, mpi_communicator);
+        for (const auto& m : rm) {
+            packReservoirConvergenceMetric(m, buf, offset, mpi_communicator);
+        }
         // Well failures.
         const auto wf = local_report.wellFailures();
         int num_wf = wf.size();
@@ -84,13 +105,16 @@ namespace
     {
         int int_pack_size = 0;
         MPI_Pack_size(1, MPI_INT, mpi_communicator, &int_pack_size);
+        int double_pack_size = 0;
+        MPI_Pack_size(1, MPI_DOUBLE, mpi_communicator, &double_pack_size);
         const int num_rf = local_report.reservoirFailures().size();
+        const int num_rm = local_report.reservoirConvergence().size();
         const int num_wf = local_report.wellFailures().size();
         int wellnames_length = 0;
         for (const auto& f : local_report.wellFailures()) {
             wellnames_length += (f.wellName().size() + 1);
         }
-        return (2 + 3*num_rf + 4*num_wf) * int_pack_size + wellnames_length;
+        return (3 + 3*num_rf + 2*num_rm + 4*num_wf)*int_pack_size + (1 + 1*num_rm)*double_pack_size + wellnames_length;
     }
 
     ConvergenceReport::ReservoirFailure unpackReservoirFailure(const std::vector<char>& recv_buffer, int& offset, MPI_Comm mpi_communicator)
@@ -107,6 +131,19 @@ namespace
                                                    phase);
     }
 
+    ConvergenceReport::ReservoirConvergenceMetric
+    unpackReservoirConvergenceMetric(const std::vector<char>& recv_buffer, int& offset, MPI_Comm mpi_communicator)
+    {
+        int type = -1;
+        int phase = -1;
+        double value = -1.0;
+        auto* data = const_cast<char*>(recv_buffer.data());
+        MPI_Unpack(data, recv_buffer.size(), &offset, &type, 1, MPI_INT, mpi_communicator);
+        MPI_Unpack(data, recv_buffer.size(), &offset, &phase, 1, MPI_INT, mpi_communicator);
+        MPI_Unpack(data, recv_buffer.size(), &offset, &value, 1, MPI_DOUBLE, mpi_communicator);
+        return { static_cast<ConvergenceReport::ReservoirFailure::Type>(type), phase, value };
+    }
+
     ConvergenceReport::WellFailure unpackWellFailure(const std::vector<char>& recv_buffer, int& offset, MPI_Comm mpi_communicator)
     {
         int type = -1;
@@ -129,14 +166,21 @@ namespace
 
     ConvergenceReport unpackSingleConvergenceReport(const std::vector<char>& recv_buffer, int& offset, MPI_Comm mpi_communicator)
     {
-        ConvergenceReport cr;
-        int num_rf = -1;
         auto* data = const_cast<char*>(recv_buffer.data());
+        double reportTime{0.0};
+        MPI_Unpack(data, recv_buffer.size(), &offset, &reportTime, 1, MPI_DOUBLE, mpi_communicator);
+        ConvergenceReport cr{reportTime};
+        int num_rf = -1;
         MPI_Unpack(data, recv_buffer.size(), &offset, &num_rf, 1, MPI_INT, mpi_communicator);
         for (int rf = 0; rf < num_rf; ++rf) {
             ConvergenceReport::ReservoirFailure f = unpackReservoirFailure(recv_buffer, offset, mpi_communicator);
             cr.setReservoirFailed(f);
         }
+        int num_rm = -1;
+        MPI_Unpack(data, recv_buffer.size(), &offset, &num_rm, 1, MPI_INT, mpi_communicator);
+        for (int rm = 0; rm < num_rm; ++rm) {
+            cr.setReservoirConvergenceMetric(unpackReservoirConvergenceMetric(recv_buffer, offset, mpi_communicator));
+        }
         int num_wf = -1;
         MPI_Unpack(data, recv_buffer.size(), &offset, &num_wf, 1, MPI_INT, mpi_communicator);
         for (int wf = 0; wf < num_wf; ++wf) {