Merge pull request #5134 from lisajulia/feature/NONNC

Feature/NONNC

CMakeLists_files.cmake

@@ -291,6 +291,10 @@ list (APPEND TEST_SOURCE_FILES
   tests/test_wellstate.cpp
   )
 
+if (HAVE_ECL_INPUT)
+  list(APPEND TEST_SOURCE_FILES tests/test_nonnc.cpp)
+endif()
+
 if(MPI_FOUND)
   list(APPEND TEST_SOURCE_FILES tests/test_parallelistlinformation.cpp
                                 tests/test_ParallelSerialization.cpp)

ebos/ecltransmissibility.hh

@@ -153,7 +153,7 @@ public:
 protected:
     void updateFromEclState_(bool global);
 
-    void removeSmallNonCartesianTransmissibilities_();
+    void removeNonCartesianTransmissibilities_(bool removeAll);
 
     /// \brief Apply the Multipliers for the case PINCH(4)==TOPBOT
     ///
@@ -287,6 +287,12 @@ protected:
     const LookUpCartesianData<Grid,GridView> lookUpCartesianData_;
 };
 
+namespace details {
+    std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2);
+    std::pair<std::uint32_t, std::uint32_t> isIdReverse(const std::uint64_t& id);
+    std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2);
+}
+
 } // namespace Opm
 
 #endif

ebos/ecltransmissibility_impl.hh

@@ -53,38 +53,37 @@
 #include <utility>
 #include <vector>
 
-namespace {
-
-constexpr unsigned elemIdxShift = 32; // bits
-
-std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
-{
-    std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
-    std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
-
-    return (elemBIdx<<elemIdxShift) + elemAIdx;
-}
-
-std::pair<std::uint32_t, std::uint32_t> isIdReverse(const std::uint64_t& id)
-{
-    // Assigning an unsigned integer to a narrower type discards the most significant bits.
-    // See "The C programming language", section A.6.2.
-    // NOTE that the ordering of element A and B may have changed
-    std::uint32_t elemAIdx = id;
-    std::uint32_t elemBIdx = (id - elemAIdx) >> elemIdxShift;
-
-    return std::make_pair(elemAIdx, elemBIdx);
-}
-
-std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
-{
-    return (std::uint64_t(elemIdx1)<<elemIdxShift) + elemIdx2;
-}
-
-}
-
 namespace Opm {
 
+namespace details {
+
+    constexpr unsigned elemIdxShift = 32; // bits
+
+    std::uint64_t isId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
+    {
+        std::uint32_t elemAIdx = std::min(elemIdx1, elemIdx2);
+        std::uint64_t elemBIdx = std::max(elemIdx1, elemIdx2);
+
+        return (elemBIdx<<elemIdxShift) + elemAIdx;
+    }
+
+    std::pair<std::uint32_t, std::uint32_t> isIdReverse(const std::uint64_t& id)
+    {
+        // Assigning an unsigned integer to a narrower type discards the most significant bits.
+        // See "The C programming language", section A.6.2.
+        // NOTE that the ordering of element A and B may have changed
+        std::uint32_t elemAIdx = id;
+        std::uint32_t elemBIdx = (id - elemAIdx) >> elemIdxShift;
+
+        return std::make_pair(elemAIdx, elemBIdx);
+    }
+
+    std::uint64_t directionalIsId(std::uint32_t elemIdx1, std::uint32_t elemIdx2)
+    {
+        return (std::uint64_t(elemIdx1)<<elemIdxShift) + elemIdx2;
+    }
+}
+
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
 EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 EclTransmissibility(const EclipseState& eclState,
@@ -114,7 +113,7 @@ template<class Grid, class GridView, class ElementMapper, class CartesianIndexMa
 Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 transmissibility(unsigned elemIdx1, unsigned elemIdx2) const
 {
-    return trans_.at(isId(elemIdx1, elemIdx2));
+    return trans_.at(details::isId(elemIdx1, elemIdx2));
 }
 
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
@@ -128,7 +127,7 @@ template<class Grid, class GridView, class ElementMapper, class CartesianIndexMa
 Scalar EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
 thermalHalfTrans(unsigned insideElemIdx, unsigned outsideElemIdx) const
 {
-    return thermalHalfTrans_.at(directionalIsId(insideElemIdx, outsideElemIdx));
+    return thermalHalfTrans_.at(details::directionalIsId(insideElemIdx, outsideElemIdx));
 }
 
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
@@ -145,7 +144,7 @@ diffusivity(unsigned elemIdx1, unsigned elemIdx2) const
     if (diffusivity_.empty())
         return 0.0;
 
-    return diffusivity_.at(isId(elemIdx1, elemIdx2));
+    return diffusivity_.at(details::isId(elemIdx1, elemIdx2));
 
 }
 
@@ -156,7 +155,7 @@ dispersivity(unsigned elemIdx1, unsigned elemIdx2) const
     if (dispersivity_.empty())
         return 0.0;
 
-    return dispersivity_.at(isId(elemIdx1, elemIdx2));
+    return dispersivity_.at(details::isId(elemIdx1, elemIdx2));
 
 }
 
@@ -175,6 +174,8 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
     const bool updateDiffusivity = eclState_.getSimulationConfig().isDiffusive();
     const bool updateDispersivity = eclState_.getSimulationConfig().rock_config().dispersion();
 
+    const bool disableNNC = eclState_.getSimulationConfig().useNONNC();
+
     if (map)
         extractPermeability_(map);
     else
@@ -334,17 +335,17 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
                 // NNC. Set zero transmissibility, as it will be
                 // *added to* by applyNncToGridTrans_() later.
                 assert(outsideFaceIdx == -1);
-                trans_[isId(elemIdx, outsideElemIdx)] = 0.0;
+                trans_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
                 if (enableEnergy_){
-                    thermalHalfTrans_[directionalIsId(elemIdx, outsideElemIdx)] = 0.0;
+                    thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = 0.0;
-                    thermalHalfTrans_[directionalIsId(outsideElemIdx, elemIdx)] = 0.0;
+                    thermalHalfTrans_[details::directionalIsId(outsideElemIdx, elemIdx)] = 0.0;
                 }
 
                 if (updateDiffusivity) {
-                    diffusivity_[isId(elemIdx, outsideElemIdx)] = 0.0;
+                    diffusivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
                 }
                 if (updateDispersivity) {
-                    dispersivity_[isId(elemIdx, outsideElemIdx)] = 0.0;
+                    dispersivity_[details::isId(elemIdx, outsideElemIdx)] = 0.0;
                 }
                 continue;
             }
@@ -457,7 +458,7 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
                                                    outsideCartElemIdx,
                                                    faceDir);
 
-            trans_[isId(elemIdx, outsideElemIdx)] = trans;
+            trans_[details::isId(elemIdx, outsideElemIdx)] = trans;
 
             // update the "thermal half transmissibility" for the intersection
             if (enableEnergy_) {
@@ -480,8 +481,8 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
                                                         axisCentroids),
                                         1.0);
                 //TODO Add support for multipliers
-                thermalHalfTrans_[directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1;
+                thermalHalfTrans_[details::directionalIsId(elemIdx, outsideElemIdx)] = halfDiffusivity1;
-                thermalHalfTrans_[directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2;
+                thermalHalfTrans_[details::directionalIsId(outsideElemIdx, elemIdx)] = halfDiffusivity2;
            }
 
             // update the "diffusive half transmissibility" for the intersection
@@ -517,7 +518,7 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
                     diffusivity = 1.0 / (1.0/halfDiffusivity1 + 1.0/halfDiffusivity2);
 
 
-                diffusivity_[isId(elemIdx, outsideElemIdx)] = diffusivity;
+                diffusivity_[details::isId(elemIdx, outsideElemIdx)] = diffusivity;
            }
 
            // update the "dispersivity half transmissibility" for the intersection
@@ -553,7 +554,7 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
                     dispersivity = 1.0 / (1.0/halfDispersivity1 + 1.0/halfDispersivity2);
 
 
-                dispersivity_[isId(elemIdx, outsideElemIdx)] = dispersivity;
+                dispersivity_[details::isId(elemIdx, outsideElemIdx)] = dispersivity;
            }
         }
     }
@@ -571,16 +572,18 @@ update(bool global, const std::function<unsigned int(unsigned int)>& map, const
         globalToLocal[cartElemIdx] = elemIdx;
     }
 
-    this->applyEditNncToGridTrans_(globalToLocal);
+    if (!disableNNC) {
-    this->applyNncToGridTrans_(globalToLocal);
+        this->applyEditNncToGridTrans_(globalToLocal);
-    this->applyEditNncrToGridTrans_(globalToLocal);
+        this->applyNncToGridTrans_(globalToLocal);
-
+        this->applyEditNncrToGridTrans_(globalToLocal);
-    if (applyNncMultregT) {
+        if (applyNncMultregT) {
-        this->applyNncMultreg_(globalToLocal);
+            this->applyNncMultreg_(globalToLocal);
+        }
     }
 
-    // Remove very small non-neighbouring transmissibilities.
+    // If disableNNC == true, remove all non-neighbouring transmissibilities.
-    this->removeSmallNonCartesianTransmissibilities_();
+    // If disableNNC == false, remove very small non-neighbouring transmissibilities.
+    this->removeNonCartesianTransmissibilities_(disableNNC);
 }
 
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
@@ -699,13 +702,14 @@ extractDispersion_()
 
 template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
 void EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>::
-removeSmallNonCartesianTransmissibilities_()
+removeNonCartesianTransmissibilities_(bool removeAll)
 {
     const auto& cartDims = cartMapper_.cartesianDimensions();
     for (auto&& trans: trans_) {
-        if (trans.second < transmissibilityThreshold_) {
+        //either remove all NNC transmissibilities or those less than the threshold (by default 1e-6 in the deck's unit system)
+        if (removeAll or trans.second < transmissibilityThreshold_) {
             const auto& id = trans.first;
-            const auto& elements = isIdReverse(id);
+            const auto& elements = details::isIdReverse(id);
             int gc1 = std::min(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
             int gc2 = std::max(cartMapper_.cartesianIndex(elements.first), cartMapper_.cartesianIndex(elements.second));
 
@@ -715,7 +719,6 @@ removeSmallNonCartesianTransmissibilities_()
             if (gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0)
                 continue;
 
-            //remove transmissibilities less than the threshold (by default 1e-6 in the deck's unit system)
             trans.second = 0.0;
         }
     }
@@ -844,7 +847,7 @@ createTransmissibilityArrays_(const std::array<bool,3>& is_tran)
             if (c1 > c2)
                 continue; // we only need to handle each connection once, thank you.
 
-            auto isID = isId(c1, c2);
+            auto isID = details::isId(c1, c2);
 
             // For CpGrid with LGRs, when leaf grid view cells with indices c1 and c2
             // have the same parent cell on level zero, then gc2 - gc1 == 0. In that case,
@@ -907,7 +910,7 @@ resetTransmissibilityFromArrays_(const std::array<bool,3>& is_tran,
             if (c1 > c2)
                 continue; // we only need to handle each connection once, thank you.
 
-            auto isID = isId(c1, c2);
+            auto isID = details::isId(c1, c2);
 
             // For CpGrid with LGRs, when leaf grid view cells with indices c1 and c2
             // have the same parent cell on level zero, then gc2 - gc1 == 0. In that case,
@@ -1012,7 +1015,7 @@ applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompr
         }
 
         {
-            auto candidate = trans_.find(isId(low, high));
+            auto candidate = trans_.find(details::isId(low, high));
             if (candidate != trans_.end()) {
                 // NNC is represented by the grid and might be a neighboring connection
                 // In this case the transmissibilty is added to the value already
@@ -1021,14 +1024,14 @@ applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompr
             }
         }
         // if (enableEnergy_) {
-        //     auto candidate = thermalHalfTrans_.find(directionalIsId(low, high));
+        //     auto candidate = thermalHalfTrans_.find(details::directionalIsId(low, high));
         //     if (candidate != trans_.end()) {
         //         // NNC is represented by the grid and might be a neighboring connection
         //         // In this case the transmissibilty is added to the value already
         //         // set or computed.
         //         candidate->second += nncEntry.transEnergy1;
         //     }
-        //     auto candidate = thermalHalfTrans_.find(directionalIsId(high, low));
+        //     auto candidate = thermalHalfTrans_.find(details::directionalIsId(high, low));
         //     if (candidate != trans_.end()) {
         //         // NNC is represented by the grid and might be a neighboring connection
         //         // In this case the transmissibilty is added to the value already
@@ -1037,7 +1040,7 @@ applyNncToGridTrans_(const std::unordered_map<std::size_t,int>& cartesianToCompr
         //     }
         // }
         // if (enableDiffusivity_) {
-        //     auto candidate = diffusivity_.find(isId(low, high));
+        //     auto candidate = diffusivity_.find(details::isId(low, high));
         //     if (candidate != trans_.end()) {
         //         // NNC is represented by the grid and might be a neighboring connection
         //         // In this case the transmissibilty is added to the value already
@@ -1127,7 +1130,7 @@ applyEditNncToGridTransHelper_(const std::unordered_map<std::size_t,int>& global
         if (low > high)
             std::swap(low, high);
 
-        auto candidate = trans_.find(isId(low, high));
+        auto candidate = trans_.find(details::isId(low, high));
         if (candidate == trans_.end()) {
             print_warning(*nnc);
             ++nnc;
@@ -1190,7 +1193,7 @@ applyNncMultreg_(const std::unordered_map<std::size_t,int>& cartesianToCompresse
                 std::swap(low, high);
             }
 
-            auto candidate = this->trans_.find(isId(low, high));
+            auto candidate = this->trans_.find(details::isId(low, high));
             if (candidate != this->trans_.end()) {
                 candidate->second *= transMult.getRegionMultiplierNNC(c1, c2);
             }

opm/simulators/utils/UnsupportedFlowKeywords.cpp

@@ -422,7 +422,6 @@ const KeywordValidation::UnsupportedKeywords& unsupportedKeywords()
         {"NOGGF", {true, std::nullopt}},
         {"NOINSPEC", {false, std::nullopt}},
         {"NOMONITO", {true, std::nullopt}},
-        {"NONNC", {true, std::nullopt}},
         {"NORSSPEC", {false, std::nullopt}},
         {"NOWARN", {false, std::nullopt}},
         {"NSTACK", {false, std::nullopt}},

tests/test_nonnc.cpp

@@ -0,0 +1,141 @@
+/*
+  Copyright 2024 Dr. Blatt - HPC-Simulation-Software & Services
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+
+#define BOOST_TEST_MODULE NONNCTest
+#define BOOST_TEST_NO_MAIN
+
+#include <opm/input/eclipse/Deck/Deck.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/EclipseGrid.hpp>
+#include <opm/input/eclipse/Parser/Parser.hpp>
+#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <opm/grid/CpGrid.hpp>
+
+#include <ebos/ecltransmissibility.hh>
+
+using namespace Opm;
+
+const int dimWorld = 3;
+const int cartDims[3] = {8,15,3};
+
+// Class extending EclTransmissibility, such that we can access the protected member trans_ to check its contents
+template<class Grid, class GridView, class ElementMapper, class CartesianIndexMapper, class Scalar>
+class Transmissibility : public EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>
+{
+    using ParentType = EclTransmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,Scalar>;
+    public:
+        Transmissibility(const EclipseState& eclState,
+                        const GridView& gridView,
+                        const CartesianIndexMapper& cartMapper,
+                        const Grid& grid,
+                        std::function<std::array<double,dimWorld>(int)> centroids,
+                        bool enableEnergy,
+                        bool enableDiffusivity,
+                        bool enableDispersivity) : ParentType(eclState,gridView,cartMapper,grid,centroids,enableEnergy,enableDiffusivity,enableDispersivity) {}
+        auto getTransmissibilitymap() {
+            return this->trans_;
+        }
+};
+
+int main(int argc, char** argv )
+{
+    Dune::MPIHelper::instance(argc, argv);
+
+    using Grid = Dune::CpGrid;
+    using GridView = Grid::LeafGridView;
+    using ElementMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GridView>;
+    using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
+    using Transmissibility = Transmissibility<Grid,GridView,ElementMapper,CartesianIndexMapper,double>;
+
+    Parser parser;
+
+    // Deck with 2 NNCs
+    const auto deck = parser.parseString( R"(RUNSPEC
+                                            NONNC
+                                            DIMENS
+                                             8 15 3 /
+                                            GRID
+                                            DX
+                                              360*10./
+                                            DY
+                                               360*10./
+                                            DZ
+                                               360*1./
+                                            TOPS
+                                               360*100./
+                                            
+                                            PORO
+                                               0. 0.25 0. 357*0.25/
+                                            PERMX
+                                                360*1000./
+                                            PERMY
+                                                360*1000./
+                                            PERMZ
+                                                360*10./
+
+                                            NNC
+                                            -- I1 J1 K1  I2 J2 K2 Trans
+                                                1  1  1   2  2  2  1000.0 / --- connection between 0 and 129
+                                                1  1  1   3  3  3  1000.0 / --- connection between 0 and 258
+                                            / 
+                                            
+                                            END)");
+    Grid grid;
+    EclipseGrid eclGrid(deck);
+    EclipseState eclState(deck);
+
+    grid.processEclipseFormat(&eclGrid, &eclState, false, false, false);
+    const auto& gridView = grid.leafGridView();
+    
+    CartesianIndexMapper cartMapper =  Dune::CartesianIndexMapper<Grid>(grid);
+
+    auto centroids = [](int) { return std::array<double,Dune::CpGrid::dimensionworld>{}; };
+    Transmissibility eclTransmissibility(eclState,
+                                         gridView,
+                                         cartMapper,
+                                         grid,
+                                         centroids,
+                                         false,false,false);
+    // Call update, true indicates that update is called on all processes
+    eclTransmissibility.update(true);
+
+    auto transmissibilityMap = eclTransmissibility.getTransmissibilitymap();
+
+    // Check if the transmissibilities of the NNCs that were added manually are either not contained in the transmissibility array (because they might be on a different process) or 0.0
+    if (transmissibilityMap.count(details::isId(0,129)) > 0)
+        assert(eclTransmissibility.transmissibility(0,129) == 0.0);
+    if (transmissibilityMap.count(details::isId(0,258)) > 0)
+        assert(eclTransmissibility.transmissibility(0,258) == 0.0);
+
+    // If there is a non-zero transmissibility in the map, ensure that it is form a neighboring connection
+    for (auto&& trans : transmissibilityMap) {
+        if (trans.second != 0.0) {
+            const auto& id = trans.first;
+            const auto& elements = details::isIdReverse(id);
+            int gc1 = std::min(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second));
+            int gc2 = std::max(cartMapper.cartesianIndex(elements.first), cartMapper.cartesianIndex(elements.second));
+            assert(gc2 - gc1 == 1 || gc2 - gc1 == cartDims[0] || gc2 - gc1 == cartDims[0]*cartDims[1] || gc2 - gc1 == 0);
+        }
+    }
+    return 0;
+}