Support for CONV_NEW when CONV in RPTRST

opm/simulators/flow/BlackoilModel.hpp

@@ -213,6 +213,14 @@ public:
     std::pair<std::vector<double>, std::vector<int>>
     characteriseCnvPvSplit(const std::vector<Scalar>& B_avg, const double dt);
 
+    /// \brief Compute the number of Newtons required by each cell in order to
+    /// satisfy the solution change convergence criteria at the last time step.
+    void convergencePerCell(const std::vector<Scalar>& B_avg,
+                            const double tol_cnv,
+                            const double dt,
+                            const double tol_cnv_energy,
+                            const int iteration);
+
     void updateTUNING(const Tuning& tuning);
 
     ConvergenceReport

opm/simulators/flow/BlackoilModel_impl.hpp

@@ -894,6 +894,9 @@ getReservoirConvergence(const double reportTime,
         }
     }
 
+    // Compute the Newton convergence per cell.
+    this->convergencePerCell(B_avg, dt, tol_cnv, tol_cnv_energy, iteration);
+
     // Output of residuals.
     if (this->terminal_output_) {
         // Only rank 0 does print to std::cout
@@ -936,6 +939,53 @@ getReservoirConvergence(const double reportTime,
     return report;
 }
 
+template <class TypeTag>
+void
+BlackoilModel<TypeTag>::
+convergencePerCell(const std::vector<Scalar>& B_avg,
+                   const double dt,
+                   const double tol_cnv,
+                   const double tol_cnv_energy,
+                   const int iteration)
+{
+    auto& rst_conv = simulator_.problem().eclWriter().mutableOutputModule().getConv();
+    if (!rst_conv.hasConv()) {
+        return;
+    }
+
+    if (iteration == 0) {
+        rst_conv.prepareConv();
+    }
+
+    const auto& residual = simulator_.model().linearizer().residual();
+    const auto& gridView = this->simulator().gridView();
+    const IsNumericalAquiferCell isNumericalAquiferCell(gridView.grid());
+    ElementContext elemCtx(this->simulator());
+    std::vector<int> convNewt(residual.size(), 0);
+    OPM_BEGIN_PARALLEL_TRY_CATCH();
+    unsigned idx = 0;
+    const int numComp = B_avg.size();
+    for (const auto& elem : elements(gridView, Dune::Partitions::interior)) {
+        elemCtx.updatePrimaryStencil(elem);
+
+        const unsigned cell_idx = elemCtx.globalSpaceIndex(/*spaceIdx=*/0, /*timeIdx=*/0);
+        const auto pvValue = simulator_.problem().referencePorosity(cell_idx, /*timeIdx=*/0) *
+                             simulator_.model().dofTotalVolume(cell_idx);
+        for (int compIdx = 0; compIdx < numComp; ++compIdx) {
+            const auto tol = (has_energy_ && compIdx == contiEnergyEqIdx) ? tol_cnv_energy : tol_cnv;
+            const Scalar cnv = std::abs(B_avg[compIdx] * residual[cell_idx][compIdx]) * dt / pvValue;
+            if (std::isnan(cnv) || cnv > maxResidualAllowed() || cnv < 0.0 || cnv > tol) {
+                convNewt[idx] = 1;
+                break;
+            }
+        }
+        ++idx;
+    }
+    OPM_END_PARALLEL_TRY_CATCH("BlackoilModel::convergencePerCell() failed: ",
+                               this->grid_.comm());
+    rst_conv.updateNewton(convNewt);
+}
+
 template <class TypeTag>
 ConvergenceReport
 BlackoilModel<TypeTag>::

opm/simulators/flow/RSTConv.cpp

@@ -42,6 +42,7 @@ void RSTConv::init(const std::size_t numCells,
     if (kw == rst_config.keywords.end()) {
         N_ = 0;
         cnv_X_.clear();
+        conv_new_.clear();
         return;
     }
 
@@ -54,6 +55,7 @@ void RSTConv::init(const std::size_t numCells,
             cnv_X_[i].resize(numCells);
         }
     }
+    conv_new_.resize(numCells, 1);
 }
 
 void RSTConv::outputRestart(data::Solution& sol)
@@ -69,6 +71,42 @@ void RSTConv::outputRestart(data::Solution& sol)
                           }
                           ++i;
                       });
+        sol.insert("CONV_NEW", std::move(conv_new_), data::TargetType::RESTART_SOLUTION);
+    }
+}
+
+bool RSTConv::hasConv() const
+{
+    return !conv_new_.empty();
+}
+
+void RSTConv::prepareConv()
+{
+    std::fill(conv_new_.begin(), conv_new_.end(), 1);
+}
+
+void RSTConv::updateNewton(const std::vector<int>& convNewt)
+{
+    const std::size_t numGloCells = conv_new_.size();
+
+    if (comm_.size() == 1) {
+        for (std::size_t n = 0; n < numGloCells; ++n) {
+            conv_new_[n] += convNewt[n];
+        }
+        return;
+    }
+
+    std::vector<int> values(numGloCells, 0);
+    const std::size_t numLocCells = convNewt.size();
+
+    for (std::size_t i = 0; i < numLocCells; ++i) {
+        values[globalCell_(i)] = convNewt[i];
+    }
+
+    comm_.sum(values.data(), numGloCells);
+
+    for (std::size_t n = 0; n < numGloCells; ++n) {
+        conv_new_[n] += values[n];
     }
 }
 

opm/simulators/flow/RSTConv.hpp

@@ -60,10 +60,23 @@ public:
     template<class ResidualVector>
     void update(const ResidualVector& resid);
 
-    //! \brief Obtain a const-ref to the accumulated data.
+    //! \brief Adds the CONV_NEW output.
+    void updateNewton(const std::vector<int>& convNewt);
+
+    //! \brief Obtain a const-ref to the accumulated data (only used for unit testing).
     const std::vector<std::vector<int>>& getData() const
     { return cnv_X_; }
 
+    //! \brief Obtain a const-ref to the CONV_NEW (only used for unit testing).
+    const std::vector<int>& getConvNew() const
+    { return conv_new_; }
+
+    //! \brief Check if CONV is required.
+    bool hasConv() const;
+
+    //! \brief CONV_NEW is restarted at iteration 0.
+    void prepareConv();
+
 private:
     //! \brief Gathers and accumulates results to the CONV arrays.
     //! \param lIdx Vector of local indices (N first is used)
@@ -77,6 +90,7 @@ private:
     LocalToGlobalCellFunc globalCell_; //!< Global cell indices
     Parallel::Communication comm_; //!< Communicator
     std::vector<std::vector<int>> cnv_X_{}; //!< Counts of worst cells for RPTRST CONV
+    std::vector<int> conv_new_{}; //!< Number of Newtons by each cell
     std::array<int,6> compIdx_{}; //!< Component indices
     int N_ = 0; //!< Number of cells to consider
 };

tests/test_rstconv.cpp

@@ -176,6 +176,17 @@ BOOST_AUTO_TEST_CASE(RstConvTest)
             }
         }
     }
+
+    std::vector<int> conv_new(10, 0);
+    cnv.updateNewton(conv_new);
+    std::vector<int> conv_new_ite1(10, 1);
+    cnv.updateNewton(conv_new_ite1);
+    if (cc.rank() == 0) {
+        BOOST_CHECK_EQUAL(cnv.getConvNew().size(), 10*cc.size());
+        for (int i = 0; i < cnv.getConvNew().size(); ++i) {
+            BOOST_CHECK_EQUAL(cnv.getConvNew()[i], 2);
+        }
+    }
 #if BOOST_VERSION / 100000 == 1 && BOOST_VERSION / 100 % 1000 < 67
     }
 #endif