From b171cbc3afe2544ae67a4e0a92b5922ac4cd0ab5 Mon Sep 17 00:00:00 2001
From: Paul <paul.egberts@tno.nl>
Date: Wed, 26 Apr 2023 17:30:53 +0200
Subject: [PATCH] Calculate headI/headJ from trajectory data when they are
 defaulted in WELSPECS

---
 .../eclipse/Schedule/Well/WellConnections.hpp |  2 ++
 .../eclipse/Schedule/KeywordHandlers.cpp      |  2 ++
 .../eclipse/Schedule/Well/WellConnections.cpp | 21 +++++++++++++++++--
 3 files changed, 23 insertions(+), 2 deletions(-)

diff --git a/opm/input/eclipse/Schedule/Well/WellConnections.hpp b/opm/input/eclipse/Schedule/Well/WellConnections.hpp
index 811478205..f61c78eec 100644
--- a/opm/input/eclipse/Schedule/Well/WellConnections.hpp
+++ b/opm/input/eclipse/Schedule/Well/WellConnections.hpp
@@ -93,6 +93,8 @@ namespace Opm {
 
         void loadWELTRAJ(const DeckRecord& record, const ScheduleGrid& grid, const std::string& wname, const KeywordLocation& location);
 
+        int getHeadI() const;
+        int getHeadJ() const;
         void add(Connection);
         std::size_t size() const;
         bool empty() const;
diff --git a/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp b/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp
index 678020436..5e5aea83e 100644
--- a/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp
+++ b/src/opm/input/eclipse/Schedule/KeywordHandlers.cpp
@@ -238,6 +238,8 @@ namespace {
                 auto connections = std::make_shared<WellConnections>(WellConnections(well2.getConnections()));
                 // cellsearchTree is calculated only once and is used to calculated cell intersections of the perforations specified in COMPTRAJ 
                 connections->loadCOMPTRAJ(record, handlerContext.grid, name, handlerContext.keyword.location(), cellSearchTree);
+                // In the case that defaults are used in WELSPECS for headI/J the headI/J are calculated based on the well trajectory data
+                well2.updateHead(connections->getHeadI(), connections->getHeadJ());
                 if (well2.updateConnections(connections, handlerContext.grid)) {
                     this->snapshots.back().wells.update( well2 );
                     wells.insert( name );
diff --git a/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp b/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp
index dde6a11f1..e0bc09b38 100644
--- a/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp
+++ b/src/opm/input/eclipse/Schedule/Well/WellConnections.cpp
@@ -577,9 +577,18 @@ namespace Opm {
             I = ijk[0];
             J = ijk[1];
             k = ijk[2];
-            // std::cout<< "I: " << I << " J: " << J << " K: " << k << std::endl;
-            external::cvf::Vec3d connection_vector = intersections[is].intersectionLengthsInCellCS;
 
+            // When using WELTRAJ & COMPTRAJ one may use default settings in WELSPECS for headI/J and let the 
+            // headI/J be calculated by the trajectory data.
+            // If these defaults are used the headI/J are set to the first intersection.
+            if (is == 0) {
+                if (this->headI == -1)
+                    this->headI = I;
+                if (this->headJ == -1)
+                    this->headJ = J;
+            }
+
+            external::cvf::Vec3d connection_vector = intersections[is].intersectionLengthsInCellCS;
 
             const CompletedCells::Cell& cell = grid.get_cell(I, J, k);
 
@@ -920,6 +929,14 @@ namespace Opm {
         return perf_length;
     }
 
+    int WellConnections::getHeadI() const {
+        return this->headI;
+    }
+
+    int WellConnections::getHeadJ() const {
+        return this->headJ;
+    }
+
     std::optional<int>
     getCompletionNumberFromGlobalConnectionIndex(const WellConnections& connections,
                                                  const std::size_t      global_index)