From b84c5f46bf5af466c2b36197748a4987f76ef3a3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= <Bard.Skaflestad@sintef.no>
Date: Mon, 28 Aug 2023 16:28:36 +0200
Subject: [PATCH] Add Method for Identifying Numerical Aquifer Cell IDs

This commit adds a new, focused, member function

  NumericalAquifers::allAquiferCellIds()

which returns a vector of those Cartesian/global cells that have
been marked as defining the model's numerical aquifers through the
AQUNUM keyword.  We intend to use this to identify those NNCs that
go to numerical aquifers--or between numerical aquifer cells--as
those may need special treatment when processing the MULTREGT
keyword.
---
 .../NumericalAquifer/NumericalAquifers.hpp    |  7 +++++--
 .../NumericalAquifer/NumericalAquifers.cpp    | 20 +++++++++++++++++++
 tests/parser/AquiferTests.cpp                 |  8 ++++++++
 3 files changed, 33 insertions(+), 2 deletions(-)

diff --git a/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.hpp b/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.hpp
index d116a8e7e..21c4fdeee 100644
--- a/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.hpp
+++ b/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.hpp
@@ -20,13 +20,15 @@
 #ifndef OPM_NUMERICALAQUIFERS_HPP
 #define OPM_NUMERICALAQUIFERS_HPP
 
+#include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
+
+#include <cstddef>
 #include <map>
 #include <unordered_map>
+#include <vector>
 
 #include <stddef.h>
 
-#include <opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/SingleNumericalAquifer.hpp>
-
 namespace Opm {
     class Deck;
     class EclipseGrid;
@@ -46,6 +48,7 @@ namespace Opm {
         bool operator==(const NumericalAquifers& other) const;
 
         std::unordered_map<size_t, const NumericalAquiferCell*> allAquiferCells() const;
+        std::vector<std::size_t> allAquiferCellIds() const;
 
         std::unordered_map<size_t, double> aquiferCellVolumes() const;
 
diff --git a/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.cpp b/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.cpp
index 29a140c8f..d8dd9f773 100644
--- a/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.cpp
+++ b/src/opm/input/eclipse/EclipseState/Aquifer/NumericalAquifer/NumericalAquifers.cpp
@@ -35,8 +35,12 @@
 
 #include <fmt/format.h>
 
+#include <algorithm>
 #include <cstddef>
 #include <unordered_set>
+#include <utility>
+#include <stdexcept>
+#include <vector>
 
 namespace Opm {
 
@@ -155,6 +159,22 @@ namespace Opm {
         return cells;
     }
 
+    std::vector<std::size_t> NumericalAquifers::allAquiferCellIds() const {
+        auto cellIds = std::vector<std::size_t>{};
+
+        for (const auto& aquifer_pair : this->m_aquifers) {
+            const auto nc = aquifer_pair.second.numCells();
+            for (auto aquCell = 0*nc; aquCell < nc; ++aquCell) {
+                cellIds.push_back(aquifer_pair.second.getCellPrt(aquCell)->global_index);
+            }
+        }
+
+        std::sort(cellIds.begin(), cellIds.end());
+        cellIds.erase(std::unique(cellIds.begin(), cellIds.end()), cellIds.end());
+
+        return cellIds;
+    }
+
     const std::map<size_t, SingleNumericalAquifer>& NumericalAquifers::aquifers() const {
         return this->m_aquifers;
     }
diff --git a/tests/parser/AquiferTests.cpp b/tests/parser/AquiferTests.cpp
index cd22fc19f..68f56510b 100644
--- a/tests/parser/AquiferTests.cpp
+++ b/tests/parser/AquiferTests.cpp
@@ -795,6 +795,14 @@ BOOST_AUTO_TEST_CASE(NumericalAquiferTest)
         BOOST_CHECK_EQUAL(c3->global_index, 3);
     }
 
+    {
+        const auto numAquCells = num_aqu.allAquiferCellIds();
+        const auto expect = std::vector { 0, 2, 3, };
+
+        BOOST_CHECK_EQUAL_COLLECTIONS(numAquCells.begin(), numAquCells.end(),
+                                      expect     .begin(), expect     .end());
+    }
+
     // using processed actnum for numerical aquifer connection generation
     std::vector<int> new_actnum(360, 1);
     new_actnum[0] = 0;