From 2964bd409e4ec141914a5621bb4ce69bfb2bd5a5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Tue, 3 Jan 2023 10:48:19 +0100
Subject: [PATCH] Use diagonal matrix block pointers for well source terms.

---
 opm/simulators/wells/BlackoilWellModel.hpp      | 2 +-
 opm/simulators/wells/BlackoilWellModel_impl.hpp | 8 ++++----
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp
index 43c426ec6..d0a761d5f 100644
--- a/opm/simulators/wells/BlackoilWellModel.hpp
+++ b/opm/simulators/wells/BlackoilWellModel.hpp
@@ -273,7 +273,7 @@ namespace Opm {
 
             // add source from wells to the reservoir matrix
             void addReservoirSourceTerms(GlobalEqVector& residual,
-                                         SparseMatrixAdapter& jacobian) const;
+                                         std::vector<typename SparseMatrixAdapter::MatrixBlock*>& diagMatAddress) const;
 
             // called at the beginning of a report step
             void beginReportStep(const int time_step);
diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index 11dbb9e85..f9daea628 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -1250,7 +1250,8 @@ namespace Opm {
 
     template <typename TypeTag>
     void BlackoilWellModel<TypeTag>::
-    addReservoirSourceTerms(GlobalEqVector& residual, SparseMatrixAdapter& jacobian) const
+    addReservoirSourceTerms(GlobalEqVector& residual,
+                            std::vector<typename SparseMatrixAdapter::MatrixBlock*>& diagMatAddress) const
     {
         // NB this loop may write multiple times to the same element
         // if a cell is perforated by more than one well, so it should
@@ -1264,14 +1265,13 @@ namespace Opm {
             for (unsigned perfIdx = 0; perfIdx < rates.size(); ++perfIdx) {
                 unsigned cellIdx = cells[perfIdx];
                 auto rate = rates[perfIdx];
-                // Scalar volume =  ebosSimulator_.problem().volume(cellIdx,0);
                 rate *= -1.0;
                 VectorBlockType res(0.0);
-		using  MatrixBlockType = Opm::MatrixBlock<Scalar, numEq, numEq >;
+                using MatrixBlockType = typename SparseMatrixAdapter::MatrixBlock;
                 MatrixBlockType bMat(0.0);
                 ebosSimulator_.model().linearizer().setResAndJacobi(res, bMat, rate);
                 residual[cellIdx] += res;
-                jacobian.addToBlock(cellIdx, cellIdx, bMat);
+                *diagMatAddress[cellIdx] += bMat;
             }
         }
     }