LevelCartesianIndexMapper for AluGrid

CMakeLists_files.cmake

@@ -780,6 +780,7 @@ list (APPEND PUBLIC_HEADER_FILES
   opm/models/utils/timerguard.hh
   opm/simulators/flow/ActionHandler.hpp
   opm/simulators/flow/AluGridCartesianIndexMapper.hpp
+  opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp
   opm/simulators/flow/AluGridVanguard.hpp
   opm/simulators/flow/Banners.hpp
   opm/simulators/flow/BaseAquiferModel.hpp

opm/simulators/flow/AluGridCartesianIndexMapper.hpp

@@ -193,10 +193,6 @@ public:
     int compressedSize() const
     { return cartesianIndex_.size(); }
 
-    /** \brief return number of cells in the active grid. Only for unifying calls with CpGrid and PolyhedralGrid specializations. */
-    int compressedLevelZeroSize() const
-    { return cartesianIndex_.size(); }
-
     /** \brief return index of the cells in the logical Cartesian grid */
     int cartesianIndex(const int compressedElementIndex) const
     {
@@ -237,15 +233,6 @@ public:
             throw std::invalid_argument("cartesianCoordinate not implemented for dimension " + std::to_string(dimension));
     }
 
-    /** \brief Only for unifying calls with CartesianIndexMapper<CpGrid> where levels are relevant */
-    void cartesianCoordinateLevel(const int compressedElementIndex, std::array<int, dimension>& coords, int level) const
-    {
-        if (level) {
-            throw std::invalid_argument("Invalid level.\n");
-        }
-        cartesianCoordinate(compressedElementIndex, coords);
-    }
-
     template <class GridView>
     std::unique_ptr<GlobalIndexDataHandle<GridView> > dataHandle(const GridView& gridView)
     {

opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp

@@ -0,0 +1,113 @@
+//===========================================================================
+//
+// File: AluGridLevelCartesianIndexMapper.hpp
+//
+// Created: Tue October 01  09:44:00 2024
+//
+// Author(s): Antonella Ritorto <antonella.ritorto@opm-op.com>
+//
+//
+// $Date$
+//
+// $Revision$
+//
+//===========================================================================
+
+/*
+  Copyright 2024 Equinor ASA.
+
+  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/>.
+*/
+#ifndef OPM_ALUGRIDLEVELCARTESIANINDEXMAPPER_HPP
+#define OPM_ALUGRIDLEVELCARTESIANINDEXMAPPER_HPP
+
+#include <dune/alugrid/grid.hh>
+#include <opm/grid/common/LevelCartesianIndexMapper.hpp>
+#include <opm/simulators/flow/AluGridCartesianIndexMapper.hpp>
+
+#include <array>
+#include <memory>
+#include <vector>
+
+
+namespace Opm
+{
+
+template<>
+class LevelCartesianIndexMapper<Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming>>
+{
+
+#if HAVE_MPI
+    using Grid = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridMPIComm>;
+#else
+    using Grid = Dune::ALUGrid<3, 3, Dune::cube, Dune::nonconforming, Dune::ALUGridNoComm>;
+#endif //HAVE_MPI
+
+ public:
+    static const int dimension = 3 ;
+
+    LevelCartesianIndexMapper(const Grid& grid,
+                              const Dune::CartesianIndexMapper<Grid>& cartesianIndexMapper)
+        : grid_{&grid}
+        , cartesianIndexMapper_{std::make_unique<Dune::CartesianIndexMapper<Grid>>(cartesianIndexMapper)}
+    {}
+          
+    const std::array<int,3>& cartesianDimensions(int level) const
+    {
+        throwIfLevelPositive(level);
+        return cartesianIndexMapper_ ->cartesianDimensions();
+    }
+
+    int cartesianSize(int level) const
+    {
+        throwIfLevelPositive(level);
+        return cartesianIndexMapper_->cartesianSize();
+    }
+
+    int compressedSize(int level) const
+    {
+        throwIfLevelPositive(level);
+        return cartesianIndexMapper_-> compressedSize();
+    }
+
+
+    int cartesianIndex( const int compressedElementIndex, const int level) const
+    {
+        throwIfLevelPositive(level);;
+        return cartesianIndexMapper_->cartesianIndex(compressedElementIndex);
+    }
+
+    void cartesianCoordinate(const int compressedElementIndex, std::array<int,dimension>& coords, int level) const
+    {
+        throwIfLevelPositive(level);
+        cartesianIndexMapper_->cartesianCoordinate(compressedElementIndex, coords);
+    }
+
+ private:
+    const Grid* grid_;
+    std::unique_ptr<Dune::CartesianIndexMapper<Grid>> cartesianIndexMapper_;
+    
+    void throwIfLevelPositive(int level) const
+    {
+        if (level) {
+            throw std::invalid_argument("Invalid level.\n");
+        }
+    }
+};
+
+}
+
+#endif

opm/simulators/flow/AluGridVanguard.hpp

@@ -34,10 +34,12 @@
 #include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <opm/grid/CpGrid.hpp>
+#include <opm/grid/cpgrid/LevelCartesianIndexMapper.hpp>
 
 #include <opm/models/common/multiphasebaseproperties.hh>
 
 #include <opm/simulators/flow/AluGridCartesianIndexMapper.hpp>
+#include <opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp>
 #include <opm/simulators/flow/FlowBaseVanguard.hpp>
 #include <opm/simulators/flow/Transmissibility.hpp>
 #include <opm/simulators/utils/ParallelEclipseState.hpp>
@@ -106,6 +108,7 @@ public:
     using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
     using GridView = GetPropType<TypeTag, Properties::GridView>;        
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
+    using LevelCartesianIndexMapper = Opm::LevelCartesianIndexMapper<Grid>;
     using EquilCartesianIndexMapper = Dune::CartesianIndexMapper<EquilGrid>;
     using TransmissibilityType = Transmissibility<Grid, GridView, ElementMapper, CartesianIndexMapper, Scalar>;
     using Factory = Dune::FromToGridFactory<Grid>;
@@ -232,6 +235,13 @@ public:
     const CartesianIndexMapper& cartesianIndexMapper() const
     { return *cartesianIndexMapper_; }
 
+    /*!
+     * \brief Returns the object which maps a global element index of the simulation grid
+     *        to the corresponding element index of the level logically Cartesian index.
+     */
+    const LevelCartesianIndexMapper levelCartesianIndexMapper() const
+    { return LevelCartesianIndexMapper(*grid_, *cartesianIndexMapper_); }
+
     /*!
      * \brief Returns mapper from compressed to cartesian indices for the EQUIL grid
      */

opm/simulators/flow/FlowProblem.hpp

@@ -253,7 +253,7 @@ public:
             // temporary measure.
             RelpermDiagnostics relpermDiagnostics{};
             relpermDiagnostics.diagnosis(simulator.vanguard().eclState(),
-                                         simulator.vanguard().cartesianIndexMapper());
+                                         simulator.vanguard().levelCartesianIndexMapper());
         }
     }
 

opm/simulators/flow/GenericCpGridVanguard.cpp

@@ -604,6 +604,13 @@ GenericCpGridVanguard<ElementMapper,GridView,Scalar>::cartesianIndexMapper() con
     return *cartesianIndexMapper_;
 }
 
+template<class ElementMapper, class GridView, class Scalar>
+const LevelCartesianIndexMapper<Dune::CpGrid>
+GenericCpGridVanguard<ElementMapper,GridView,Scalar>::levelCartesianIndexMapper() const
+{
+    return LevelCartesianIndexMapper(*grid_);
+}
+
 template<class ElementMapper, class GridView, class Scalar>
 const Dune::CartesianIndexMapper<Dune::CpGrid>&
 GenericCpGridVanguard<ElementMapper,GridView,Scalar>::equilCartesianIndexMapper() const

opm/simulators/flow/GenericCpGridVanguard.hpp

@@ -28,6 +28,7 @@
 #define OPM_GENERIC_CPGRID_VANGUARD_HPP
 
 #include <opm/grid/CpGrid.hpp>
+#include <opm/grid/cpgrid/LevelCartesianIndexMapper.hpp>
 
 #include <opm/simulators/flow/FlowGenericVanguard.hpp>
 
@@ -78,6 +79,7 @@ template<class ElementMapper, class GridView, class Scalar>
 class GenericCpGridVanguard {
 protected:
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Dune::CpGrid>;
+    using LevelCartesianIndexMapper = Opm::LevelCartesianIndexMapper<Dune::CpGrid>;
     using Element = typename GridView::template Codim<0>::Entity;
 
 public:
@@ -131,6 +133,12 @@ public:
      */
     const CartesianIndexMapper& cartesianIndexMapper() const;
 
+    /*!
+     * \brief Returns the object which maps a global element index of the simulation grid
+     *        to the corresponding element index of the level logically Cartesian index.
+     */
+    const LevelCartesianIndexMapper levelCartesianIndexMapper() const;
+
     /*!
      * \brief Returns mapper from compressed to cartesian indices for the EQUIL grid
      */

opm/simulators/flow/PolyhedralGridVanguard.hpp

@@ -28,6 +28,7 @@
 #define OPM_POLYHEDRAL_GRID_VANGUARD_HPP
 
 #include <opm/grid/polyhedralgrid.hh>
+#include <opm/grid/polyhedralgrid/levelcartesianindexmapper.hh>
 
 #include <opm/models/common/multiphasebaseproperties.hh>
 
@@ -93,6 +94,7 @@ public:
     using EquilGrid = GetPropType<TypeTag, Properties::EquilGrid>;
     using GridView = GetPropType<TypeTag, Properties::GridView>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
+    using LevelCartesianIndexMapper = Opm::LevelCartesianIndexMapper<Grid>;
     using EquilCartesianIndexMapper = Dune::CartesianIndexMapper<EquilGrid>;
     static constexpr int dimension = Grid::dimension;
     static constexpr int dimensionworld = Grid::dimensionworld;
@@ -169,6 +171,13 @@ public:
     const CartesianIndexMapper& cartesianIndexMapper() const
     { return *cartesianIndexMapper_; }
 
+    /*!
+     * \brief Returns the object which maps a global element index of the simulation grid
+     *        to the corresponding element index of the level logically Cartesian index.
+     */
+    const LevelCartesianIndexMapper levelCartesianIndexMapper() const
+    { return LevelCartesianIndexMapper(*grid_); }
+
     /*!
      * \brief Returns mapper from compressed to cartesian indices for the EQUIL grid
      *

opm/simulators/utils/satfunc/RelpermDiagnostics.cpp

@@ -45,10 +45,13 @@
 
 #include <opm/grid/CpGrid.hpp>
 #include <opm/grid/polyhedralgrid.hh>
+#include <opm/grid/cpgrid/LevelCartesianIndexMapper.hpp>
+#include <opm/grid/polyhedralgrid/levelcartesianindexmapper.hh>
 
 #ifdef HAVE_DUNE_ALUGRID
 #include <dune/alugrid/grid.hh>
 #include <dune/alugrid/3d/gridview.hh>
+#include <opm/simulators/flow/AluGridLevelCartesianIndexMapper.hpp>
 #endif // HAVE_DUNE_ALUGRID
 
 namespace Opm {
@@ -822,9 +825,9 @@ namespace Opm {
         }
     }
 
-    template <class CartesianIndexMapper>
+    template <class LevelCartesianIndexMapper>
     void RelpermDiagnostics::diagnosis(const EclipseState& eclState,
-                                       const CartesianIndexMapper& cartesianIndexMapper)
+                                       const LevelCartesianIndexMapper& levelCartesianIndexMapper)
     {
         OpmLog::info("\n===============Saturation Functions Diagnostics===============\n");
         bool doDiagnostics = phaseCheck_(eclState);
@@ -833,16 +836,16 @@ namespace Opm {
         satFamilyCheck_(eclState);
         tableCheck_(eclState);
         unscaledEndPointsCheck_(eclState);
-        scaledEndPointsCheck_(eclState, cartesianIndexMapper);
+        scaledEndPointsCheck_(eclState, levelCartesianIndexMapper);
     }
 
-    template <class CartesianIndexMapper>
+    template <class LevelCartesianIndexMapper>
     void RelpermDiagnostics::scaledEndPointsCheck_(const EclipseState& eclState,
-                                                   const CartesianIndexMapper& cartesianIndexMapper)
+                                                   const LevelCartesianIndexMapper& levelCartesianIndexMapper)
     {
         // All end points are subject to round-off errors, checks should account for it
         const float tolerance = 1e-6;
-        const int nc = cartesianIndexMapper.compressedLevelZeroSize();
+        const int nc = levelCartesianIndexMapper.compressedSize(0);
         const bool threepoint = eclState.runspec().endpointScaling().threepoint();
         scaledEpsInfo_.resize(nc);
         EclEpsGridProperties epsGridProperties(eclState, false);
@@ -852,7 +855,7 @@ namespace Opm {
             std::string cellIdx;
             {
                 std::array<int, 3> ijk;
-                cartesianIndexMapper.cartesianCoordinateLevel(c, ijk, 0);
+                levelCartesianIndexMapper.cartesianCoordinate(c, ijk, 0);
                 cellIdx = "(" + std::to_string(ijk[0]) + ", " +
                     std::to_string(ijk[1]) + ", " +
                     std::to_string(ijk[2]) + ")";
@@ -887,8 +890,8 @@ namespace Opm {
     }
 
 #define INSTANCE_DIAGNOSIS(...) \
-    template void RelpermDiagnostics::diagnosis<Dune::CartesianIndexMapper<__VA_ARGS__>>(const EclipseState&, const Dune::CartesianIndexMapper<__VA_ARGS__>&); \
+    template void RelpermDiagnostics::diagnosis<Opm::LevelCartesianIndexMapper<__VA_ARGS__>>(const EclipseState&, const Opm::LevelCartesianIndexMapper<__VA_ARGS__>&); \
-    template void RelpermDiagnostics::scaledEndPointsCheck_<Dune::CartesianIndexMapper<__VA_ARGS__>>(const EclipseState&, const Dune::CartesianIndexMapper<__VA_ARGS__>&);
+    template void RelpermDiagnostics::scaledEndPointsCheck_<Opm::LevelCartesianIndexMapper<__VA_ARGS__>>(const EclipseState&, const Opm::LevelCartesianIndexMapper<__VA_ARGS__>&);
 
     INSTANCE_DIAGNOSIS(Dune::CpGrid)
     INSTANCE_DIAGNOSIS(Dune::PolyhedralGrid<3,3>)

opm/simulators/utils/satfunc/RelpermDiagnostics.hpp

@@ -28,6 +28,7 @@
 namespace Opm {
 
     class EclipseState;
+    template <typename Grid> class LevelCartesianIndexMapper;
     class MiscTable;
     class MsfnTable;
     class SgcwmisTable;
@@ -54,9 +55,9 @@ namespace Opm {
         ///output if they're found.
         ///\param[in] eclState  eclipse state.
         ///\param[in] grid      unstructured grid.
-        template <class CartesianIndexMapper>
+        template <class LevelCartesianIndexMapper>
         void diagnosis(const EclipseState& eclState,
-                       const CartesianIndexMapper& cartesianIndexMapper);
+                       const LevelCartesianIndexMapper& levelCartesianIndexMapper);
 
     private:
         enum FluidSystem {
@@ -95,9 +96,9 @@ namespace Opm {
         ///Check endpoints in the saturation tables.
         void unscaledEndPointsCheck_(const EclipseState& eclState);
 
-        template <class CartesianIndexMapper>
+        template <class LevelCartesianIndexMapper>
         void scaledEndPointsCheck_(const EclipseState& eclState,
-                                   const CartesianIndexMapper& cartesianIndexMapper);
+                                   const LevelCartesianIndexMapper& levelCartesianIndexMapper);
 
         ///For every table, need to deal with case by case.
         void swofTableCheck_(const SwofTable& swofTables,

tests/test_relpermdiagnostics.cpp

@@ -41,6 +41,7 @@
 #include <opm/common/OpmLog/CounterLog.hpp>
 
 #include <opm/grid/CpGrid.hpp>
+#include <opm/grid/cpgrid/LevelCartesianIndexMapper.hpp>
 #include <dune/grid/common/mcmgmapper.hh>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
@@ -84,12 +85,12 @@ BOOST_AUTO_TEST_CASE(diagnosis)
                                /*flipNormals=*/false,
                                /*clipZ=*/false);
 
-    typedef Dune::CartesianIndexMapper<Grid> CartesianIndexMapper;
+    typedef Opm::LevelCartesianIndexMapper<Grid> LevelCartesianIndexMapper;
-    CartesianIndexMapper cartesianIndexMapper = CartesianIndexMapper(grid);
+    LevelCartesianIndexMapper levelCartesianIndexMapper = LevelCartesianIndexMapper(grid);
     std::shared_ptr<CounterLog> counterLog = std::make_shared<CounterLog>(Log::DefaultMessageTypes);
     OpmLog::addBackend( "COUNTERLOG" , counterLog );
     RelpermDiagnostics diagnostics;
-    diagnostics.diagnosis(eclState, cartesianIndexMapper);
+    diagnostics.diagnosis(eclState, levelCartesianIndexMapper);
     BOOST_CHECK_EQUAL(1, counterLog->numMessages(Log::MessageType::Warning));
 }
 BOOST_AUTO_TEST_SUITE_END()