dunify the aquifer models

opm/simulators/aquifers/AquiferCarterTracy.hpp

@@ -46,6 +46,7 @@ public:
     using typename Base::RateVector;
     using typename Base::Scalar;
     using typename Base::Simulator;
+    using typename Base::ElementMapper;
 
     using Base::waterCompIdx;
     using Base::waterPhaseIdx;
@@ -92,7 +93,6 @@ protected:
     inline void initializeConnections() override
     {
         const auto& eclState = this->ebos_simulator_.vanguard().eclState();
-        const auto& ugrid = this->ebos_simulator_.vanguard().grid();
         const auto& grid = eclState.getInputGrid();
 
         // We hack the cell depth values for now. We can actually get it from elementcontext pos
@@ -100,9 +100,6 @@ protected:
         this->alphai_.resize(this->size(), 1.0);
         this->faceArea_connected_.resize(this->size(), 0.0);
 
-        auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
-        auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
-
         // Translate the C face tag into the enum used by opm-parser's TransMult class
         Opm::FaceDir::DirEnum faceDirection;
 
@@ -110,18 +107,38 @@ protected:
         Scalar denom_face_areas = 0.;
         this->cellToConnectionIdx_.resize(this->ebos_simulator_.gridView().size(/*codim=*/0), -1);
         for (size_t idx = 0; idx < this->size(); ++idx) {
-            const int cell_index = this->cartesian_to_compressed_.at(this->connections_[idx].global_index);
+            const auto global_index = this->connections_[idx].global_index;
+            const int cell_index = this->cartesian_to_compressed_.at(global_index);
             this->cellToConnectionIdx_[cell_index] = idx;
+            const auto cellCenter = grid.getCellCenter(global_index);
+            this->cell_depth_.at(idx) = cellCenter[2];
+        }
+        // get default areas for all intersections
+        const auto& gridView = this->ebos_simulator_.vanguard().gridView();
+        ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
+        auto elemIt = gridView.template begin</*codim=*/ 0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const auto& elem = *elemIt;
+            unsigned cell_index = elemMapper.index(elem);
+            int idx = this->cellToConnectionIdx_[cell_index];
 
-            const auto cellFacesRange = cell2Faces[cell_index];
+            // only deal with connections given by the aquifer
-            for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
+            if( idx < 0)
-                // The index of the face in the compressed grid
+                continue;
-                const int faceIdx = *cellFaceIter;
 
-                // the logically-Cartesian direction of the face
+            auto isIt = gridView.ibegin(elem);
-                const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
+            const auto& isEndIt = gridView.iend(elem);
+            for (; isIt != isEndIt; ++ isIt) {
+                // store intersection, this might be costly
+                const auto& intersection = *isIt;
 
-                switch (faceTag) {
+                // only deal with grid boundaries
+                if (!intersection.boundary())
+                    continue;
+
+                int insideFaceIdx  = intersection.indexInInside();
+                switch (insideFaceIdx) {
                 case 0:
                     faceDirection = Opm::FaceDir::XMinus;
                     break;
@@ -146,12 +163,10 @@ protected:
                 }
 
                 if (faceDirection == this->connections_[idx].face_dir) {
-                    this->faceArea_connected_.at(idx) = this->getFaceArea(faceCells, ugrid, faceIdx, idx);
+                    this->faceArea_connected_[idx] = this->getFaceArea(intersection, idx);
                     denom_face_areas += (this->connections_[idx].influx_mult * this->faceArea_connected_.at(idx));
                 }
             }
-            auto cellCenter = grid.getCellCenter(this->connections_[idx].global_index);
-            this->cell_depth_.at(idx) = cellCenter[2];
         }
 
         const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());

opm/simulators/aquifers/AquiferFetkovich.hpp

@@ -47,6 +47,7 @@ public:
     using typename Base::RateVector;
     using typename Base::Scalar;
     using typename Base::Simulator;
+    using typename Base::ElementMapper;
 
     using Base::waterCompIdx;
     using Base::waterPhaseIdx;
@@ -95,7 +96,6 @@ protected:
     inline void initializeConnections() override
     {
         const auto& eclState = this->ebos_simulator_.vanguard().eclState();
-        const auto& ugrid = this->ebos_simulator_.vanguard().grid();
         const auto& grid = eclState.getInputGrid();
 
         // We hack the cell depth values for now. We can actually get it from elementcontext pos
@@ -103,9 +103,6 @@ protected:
         this->alphai_.resize(this->size(), 1.0);
         this->faceArea_connected_.resize(this->size(), 0.0);
 
-        auto cell2Faces = Opm::UgGridHelpers::cell2Faces(ugrid);
-        auto faceCells = Opm::UgGridHelpers::faceCells(ugrid);
-
         // Translate the C face tag into the enum used by opm-parser's TransMult class
         Opm::FaceDir::DirEnum faceDirection;
 
@@ -115,21 +112,37 @@ protected:
         for (size_t idx = 0; idx < this->size(); ++idx) {
             const auto global_index = this->connections_[idx].global_index;
             const int cell_index = this->cartesian_to_compressed_.at(global_index);
-
             this->cellToConnectionIdx_[cell_index] = idx;
             const auto cellCenter = grid.getCellCenter(global_index);
             this->cell_depth_.at(idx) = cellCenter[2];
+        }
+        // get areas for all connections
+        const auto& gridView = this->ebos_simulator_.vanguard().gridView();
+        ElementMapper elemMapper(gridView, Dune::mcmgElementLayout());
+        auto elemIt = gridView.template begin</*codim=*/ 0>();
+        const auto& elemEndIt = gridView.template end</*codim=*/ 0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const auto& elem = *elemIt;
+            unsigned cell_index = elemMapper.index(elem);
+            int idx = this->cellToConnectionIdx_[cell_index];
+
+            // only deal with connections given by the aquifer
+            if( idx < 0)
+                continue;
 
             if (!this->connections_[idx].influx_coeff.first) { // influx_coeff is defaulted
-                const auto cellFacesRange = cell2Faces[cell_index];
+                auto isIt = gridView.ibegin(elem);
-                for (auto cellFaceIter = cellFacesRange.begin(); cellFaceIter != cellFacesRange.end(); ++cellFaceIter) {
+                const auto& isEndIt = gridView.iend(elem);
-                    // The index of the face in the compressed grid
+                for (; isIt != isEndIt; ++ isIt) {
-                    const int faceIdx = *cellFaceIter;
+                    // store intersection, this might be costly
+                    const auto& intersection = *isIt;
 
-                    // the logically-Cartesian direction of the face
+                    // only deal with grid boundaries
-                    const int faceTag = Opm::UgGridHelpers::faceTag(ugrid, cellFaceIter);
+                    if (!intersection.boundary())
+                        continue;
 
-                    switch (faceTag) {
+                    int insideFaceIdx  = intersection.indexInInside();
+                    switch (insideFaceIdx) {
                     case 0:
                         faceDirection = Opm::FaceDir::XMinus;
                         break;
@@ -150,11 +163,11 @@ protected:
                         break;
                     default:
                         OPM_THROW(Opm::NumericalIssue,
-                                  "Initialization of Aquifer problem. Make sure faceTag is correctly defined");
+                            "Initialization of Aquifer Fetkovich problem. Make sure faceTag is correctly defined");                }
-                    }
+
 
                     if (faceDirection == this->connections_[idx].face_dir) {
-                        this->faceArea_connected_[idx] = this->getFaceArea(faceCells, ugrid, faceIdx, idx);
+                        this->faceArea_connected_[idx] = this->getFaceArea(intersection, idx);
                         break;
                     }
                 }

opm/simulators/aquifers/AquiferInterface.hpp

@@ -50,6 +50,7 @@ public:
     using BlackoilIndices = GetPropType<TypeTag, Properties::Indices>;
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
 
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
@@ -201,29 +202,13 @@ protected:
         rhow_.at(idx) = fs.density(waterPhaseIdx);
     }
 
-    template <class faceCellType, class ugridType>
+    template <class Intersection>
-    inline double getFaceArea(const faceCellType& faceCells,
+    inline double getFaceArea(const Intersection& intersection,
-                              const ugridType& ugrid,
+                              unsigned idx) const
-                              const int faceIdx,
-                              const int idx) const
     {
-        // Check now if the face is outside of the reservoir, or if it adjoins an inactive cell
+        const auto& geometry = intersection.geometry();
-        // Do not make the connection if the product of the two cellIdx > 0. This is because the
+        const auto defaultFaceArea = geometry.volume();
-        // face is within the reservoir/not connected to boundary. (We still have yet to check for inactive cell
+        return (!this->connections_[idx].influx_coeff.first) ? defaultFaceArea : this->connections_[idx].influx_coeff.second;
-        // adjoining)
-        double faceArea = 0.;
-        const auto cellNeighbour0 = faceCells(faceIdx, 0);
-        const auto cellNeighbour1 = faceCells(faceIdx, 1);
-        const auto defaultFaceArea = Opm::UgGridHelpers::faceArea(ugrid, faceIdx);
-        const auto calculatedFaceArea
-            = (!this->connections_[idx].influx_coeff.first) ? defaultFaceArea : this->connections_[idx].influx_coeff.second;
-        faceArea = (cellNeighbour0 * cellNeighbour1 > 0) ? 0. : calculatedFaceArea;
-        if (cellNeighbour1 == 0) {
-            faceArea = (cellNeighbour0 < 0) ? faceArea : 0.;
-        } else if (cellNeighbour0 == 0) {
-            faceArea = (cellNeighbour1 < 0) ? faceArea : 0.;
-        }
-        return faceArea;
     }
 
     virtual void endTimeStep() = 0;

opm/simulators/aquifers/BlackoilAquiferModel_impl.hpp

@@ -170,11 +170,10 @@ BlackoilAquiferModel<TypeTag>::init()
         throw std::runtime_error("Aquifers currently do not work in parallel.");
 
     // Get all the carter tracy aquifer properties data and put it in aquifers vector
-    const auto& ugrid = simulator_.vanguard().grid();
     const int number_of_cells = simulator_.gridView().size(0);
 
     cartesian_to_compressed_ = cartesianToCompressed(number_of_cells,
-                                                     Opm::UgGridHelpers::globalCell(ugrid));
+                                                     this->simulator_.vanguard().grid().globalCell().data());
 
     const auto& connections = aquifer.connections();
     for (const auto& aq : aquifer.ct()) {