using the depth specified by numerical aquifer

for numerical aquifer cells, instead of the depth calculated from the
grid geometry.

ebos/eclbasevanguard.hh

@@ -791,6 +791,15 @@ protected:
             const Element& element = *elemIt;
             const unsigned int elemIdx = elemMapper.index(element);
             cellCenterDepth_[elemIdx] = cellCenterDepth(element);
+
+            if (this->eclState_->aquifer().hasNumericalAquifer()) {
+                const auto& num_aquifer = this->eclState_->aquifer().numericalAquifers();
+                const unsigned int global_index = cartesianIndex(elemIdx);
+                if (num_aquifer.hasCell(global_index)) {
+                    const auto& cell = num_aquifer.aquiferCells().at(global_index);
+                    cellCenterDepth_[elemIdx] = cell.depth;
+                }
+            }
         }
     }
 

ebos/eclproblem.hh

@@ -2735,6 +2735,13 @@ private:
             auto& elemFluidState = initialFluidStates_[elemIdx];
             elemFluidState.assign(equilInitializer.initialFluidState(elemIdx));
         }
+        const auto& aquifer_cells = simulator.vanguard().eclState().aquifer().numericalAquifers().aquiferCells();
+
+        const auto& vanguard = simulator.vanguard();
+        size_t numCartDof = vanguard.cartesianSize();
+        std::vector<int> cartesianToCompressedElemIdx(numCartDof, -1);
+        for (unsigned elemIdx = 0; elemIdx < numElems; ++elemIdx)
+            cartesianToCompressedElemIdx[vanguard.cartesianIndex(elemIdx)] = elemIdx;
     }
 
     void readEclRestartSolution_()

ebos/equil/initstateequil.hh

@@ -1591,6 +1591,7 @@ public:
         }
 
         //Querry cell depth, cell top-bottom.
+        // aquifer should enter here
         updateCellProps_(gridView);
 
         // Get the equilibration records.
@@ -1914,7 +1915,7 @@ private:
         using CellID  = std::remove_cv_t<std::remove_reference_t<
             decltype(std::declval<CellPos>().cell)>>;
 
-        this->cellLoop(cells, [this, &eqreg, &ptable, &psat]
+        this->cellLoop(cells, [this, &eqreg,  &ptable, &psat]
             (const CellID                 cell,
              Details::PhaseQuantityValue& pressures,
              Details::PhaseQuantityValue& saturations,
@@ -1925,6 +1926,16 @@ private:
                 cell, cellCenterDepth_[cell]
             };
 
+/*            const size_t global_index = UgGridHelpers::globalCell(grid)[cell];
+
+            if (aquifer.hasCell(global_index)) {
+                const auto& aqu_cells = aquifer.aquiferCells();
+                const auto& aqu_cell = aqu_cells.at(global_index);
+                pos = CellPos {
+                    cell, aqu_cell.depth
+                };
+            } */
+
             saturations = psat.deriveSaturations(pos, eqreg, ptable);
             pressures   = psat.correctedPhasePressures();