From 91a54645fb1aa937a3336d4886ca89549ffbcedf Mon Sep 17 00:00:00 2001
From: Tong Dong Qiu <tongdong.qiu@bigdataccelerate.com>
Date: Tue, 7 Jun 2022 10:21:38 +0200
Subject: [PATCH] Fix: zeros were not replaced correctly in the diagonal And
 minor changes

---
 opm/simulators/linalg/bda/BdaBridge.cpp | 15 +++++++--------
 opm/simulators/linalg/bda/BdaBridge.hpp |  2 ++
 2 files changed, 9 insertions(+), 8 deletions(-)

diff --git a/opm/simulators/linalg/bda/BdaBridge.cpp b/opm/simulators/linalg/bda/BdaBridge.cpp
index 822561773..bd9e1d72b 100644
--- a/opm/simulators/linalg/bda/BdaBridge.cpp
+++ b/opm/simulators/linalg/bda/BdaBridge.cpp
@@ -128,17 +128,16 @@ BdaBridge<BridgeMatrix, BridgeVector, block_size>::BdaBridge(std::string acceler
 
 
 template <class BridgeMatrix>
-int checkZeroDiagonal(BridgeMatrix& mat) {
-    static std::vector<typename BridgeMatrix::size_type> diag_indices;   // contains offsets of the diagonal nnzs
+int replaceZeroDiagonal(BridgeMatrix& mat, std::vector<typename BridgeMatrix::size_type>& diag_indices) {
     int numZeros = 0;
-    const int dim = 3;                    // might be replaced with mat[0][0].N() or BridgeMatrix::block_type::size()
+    const int dim = mat[0][0].N();                    // might be replaced with BridgeMatrix::block_type::size()
     const double zero_replace = 1e-15;
     if (diag_indices.empty()) {
-        int N = mat.N();
-        diag_indices.reserve(N);
+        int Nb = mat.N();
+        diag_indices.reserve(Nb);
         for (typename BridgeMatrix::iterator r = mat.begin(); r != mat.end(); ++r) {
             auto diag = r->find(r.index());  // diag is an iterator
-            assert(diag.index() == r.index());
+            assert(diag.index() == r.index()); // every row must have a diagonal block
             for (int rr = 0; rr < dim; ++rr) {
                 auto& val = (*diag)[rr][rr]; // reference to easily change the value
                 if (val == 0.0) {             // could be replaced by '< 1e-30' or similar
@@ -244,7 +243,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system(BridgeMatri
         }
 
         Dune::Timer t_zeros;
-        int numZeros = checkZeroDiagonal(*bridgeMat);
+        int numZeros = replaceZeroDiagonal(*bridgeMat, diagIndices);
         if (verbosity >= 2) {
             std::ostringstream out;
             out << "Checking zeros took: " << t_zeros.stop() << " s, found " << numZeros << " zeros";
@@ -263,7 +262,7 @@ void BdaBridge<BridgeMatrix, BridgeVector, block_size>::solve_system(BridgeMatri
             }
 
             Dune::Timer t_zeros2;
-            int jacNumZeros = checkZeroDiagonal(*jacMat);
+            int jacNumZeros = replaceZeroDiagonal(*jacMat, jacDiagIndices);
             if (verbosity >= 2) {
                 std::ostringstream out;
                 out << "Checking zeros for jacMat took: " << t_zeros2.stop() << " s, found " << jacNumZeros << " zeros";
diff --git a/opm/simulators/linalg/bda/BdaBridge.hpp b/opm/simulators/linalg/bda/BdaBridge.hpp
index 4e79b9eae..9133c5a59 100644
--- a/opm/simulators/linalg/bda/BdaBridge.hpp
+++ b/opm/simulators/linalg/bda/BdaBridge.hpp
@@ -48,6 +48,8 @@ private:
     std::shared_ptr<Opm::Accelerator::BlockedMatrix> jacMatrix;  // 'stores' preconditioner matrix, actually points to h_rows, h_cols and the received BridgeMatrix for the nonzeroes
     std::vector<int> h_rows, h_cols;  // store the sparsity pattern of the matrix
     std::vector<int> h_jacRows, h_jacCols;  // store the sparsity pattern of the jacMatrix
+    std::vector<typename BridgeMatrix::size_type> diagIndices;   // contains offsets of the diagonal blocks wrt start of the row, used for replaceZeroDiagonal()
+    std::vector<typename BridgeMatrix::size_type> jacDiagIndices;   // same but for jacMatrix
 
 public:
     /// Construct a BdaBridge