Periodic index spaces

example/grid_tutorial/06_local_grid/local_grid_example.cpp

@@ -223,6 +223,23 @@ void localGridExample()
               << boundary_cell_space.size() << "\n"
               << std::endl;
 
+    /*
+      The directions which are shared along periodic edges of the simulation
+      domain constitute the periodic index spaces. These are the complement of
+      the boundary index space (periodic vs open system boundaries) and are a
+      subset of the shared index spaces (any periodic boundary is shared with
+      another rank, even if it is a self neighbor).
+    */
+    auto ghost_periodic_cell_space = local_grid->periodicIndexSpace(
+        Cabana::Grid::Ghost(), Cabana::Grid::Cell(), -1, 0, 1 );
+    std::cout << "Periodic index space (Ghost, Cell):\nMin: ";
+    for ( int d = 0; d < 3; ++d )
+        std::cout << ghost_periodic_cell_space.min( d ) << " ";
+    std::cout << "\nMax: ";
+    for ( int d = 0; d < 3; ++d )
+        std::cout << ghost_periodic_cell_space.max( d ) << " ";
+    std::cout << "\n" << std::endl;
+
     /*
       We now create a partially non-periodic global grid to highlight some
       details of the local grid boundary index spaces.

grid/src/Cabana_Grid_LocalGrid.hpp

@@ -251,8 +251,82 @@ class LocalGrid
     boundaryIndexSpace( DecompositionTag t1, EntityType t2, const int off_i,
                         const int off_j, const int halo_width = -1 ) const;
 
+    /*!
+      \brief Given the relative offsets of a periodic boundary relative to this
+      local grid's indices get the set of local entity indices associated with
+      that boundary in the given decomposition.
+
+      \param t1 Decomposition type: Own or Ghost
+      \param t2 Entity: Cell, Node, Edge, or Face
+      \param off_ijk %Array of neighbor offset indices.
+      \param halo_width Optional depth of shared indices within the halo. Must
+      be less than or equal to the halo width of the local grid. Default is to
+      use the halo width of the local grid.
+
+      For example, if the Own decomposition is used, the interior entities that
+      would be affected by a periodic boundary operation are provided whereas if
+      the Ghost decomposition is used the halo entities of that periodic
+      neighbor are provided.
+    */
+    template <class DecompositionTag, class EntityType>
+    IndexSpace<num_space_dim>
+    periodicIndexSpace( DecompositionTag t1, EntityType t2,
+                        const std::array<int, num_space_dim>& off_ijk,
+                        const int halo_width = -1 ) const;
+
+    /*!
+      \brief Given the relative offsets of a periodic boundary relative to this
+      local grid's indices get the set of local entity indices associated with
+      that boundary in the given decomposition.
+
+      \param t1 Decomposition type: Own or Ghost
+      \param t2 Entity: Cell, Node, Edge, or Face
+      \param off_i, off_j, off_k Neighbor offset index in a given dimension.
+      \param halo_width Optional depth of shared indices within the halo. Must
+      be less than or equal to the halo width of the local grid. Default is to
+      use the halo width of the local grid.
+
+      For example, if the Own decomposition is used, the interior entities that
+      would be affected by a boundary operation are provided whereas if the
+      Ghost decomposition is used the halo entities of that periodic neighbor
+      are provided.
+    */
+    template <class DecompositionTag, class EntityType,
+              std::size_t NSD = num_space_dim>
+    std::enable_if_t<3 == NSD, IndexSpace<3>>
+    periodicIndexSpace( DecompositionTag t1, EntityType t2, const int off_i,
+                        const int off_j, const int off_k,
+                        const int halo_width = -1 ) const;
+
+    /*!
+      \brief Given the relative offsets of a periodic boundary relative to this
+      local grid's indices get the set of local entity indices associated with
+      that boundary in the given decomposition.
+
+      \param t1 Decomposition type: Own or Ghost
+      \param t2 Entity: Cell, Node, Edge, or Face
+      \param off_i, off_j Neighbor offset index in a given dimension.
+      \param halo_width Optional depth of shared indices within the halo. Must
+      be less than or equal to the halo width of the local grid. Default is to
+      use the halo width of the local grid.
+
+      For example, if the Own decomposition is used, the interior entities that
+      would be affected by a boundary operation are provided whereas if the
+      Ghost decomposition is used the halo entities of that periodic neighbor
+      are provided.
+    */
+    template <class DecompositionTag, class EntityType,
+              std::size_t NSD = num_space_dim>
+    std::enable_if_t<2 == NSD, IndexSpace<2>>
+    periodicIndexSpace( DecompositionTag t1, EntityType t2, const int off_i,
+                        const int off_j, const int halo_width = -1 ) const;
+
   private:
     // Helper functions
+    auto zeroIndexSpace() const;
+    auto setupHaloWidth( const int halo_width ) const;
+    void checkOffsets( const std::array<int, num_space_dim>& off_ijk ) const;
+
     template <class OwnedIndexSpace>
     auto getBound( OwnedIndexSpace owned_space, const int upper_lower,
                    const std::array<int, num_space_dim>& off_ijk,

grid/src/Cabana_Grid_LocalGrid_impl.hpp

@@ -92,6 +92,42 @@ LocalGrid<MeshType>::neighborRank( const int off_i, const int off_j ) const
 }
 
 //! \cond Impl
+//---------------------------------------------------------------------------//
+// Get the index space for a given combination of decomposition, entity, and
+// index types.
+template <class MeshType>
+auto LocalGrid<MeshType>::zeroIndexSpace() const
+{
+    std::array<long, num_space_dim> zero_size;
+    for ( std::size_t d = 0; d < num_space_dim; ++d )
+        zero_size[d] = 0;
+    return IndexSpace<num_space_dim>( zero_size, zero_size );
+}
+
+template <class MeshType>
+auto LocalGrid<MeshType>::setupHaloWidth( const int halo_width ) const
+{
+    // If we got the default halo width of -1 this means we want to use the
+    // default of the entire halo.
+    int hw = ( -1 == halo_width ) ? _halo_cell_width : halo_width;
+
+    // Check that the requested halo width is valid.
+    if ( hw > _halo_cell_width )
+        throw std::logic_error(
+            "Requested halo width larger than local grid halo" );
+    return hw;
+}
+
+template <class MeshType>
+void LocalGrid<MeshType>::checkOffsets(
+    const std::array<int, num_space_dim>& off_ijk ) const
+{
+    // Check that the offsets are valid.
+    for ( std::size_t d = 0; d < num_space_dim; ++d )
+        if ( off_ijk[d] < -1 || 1 < off_ijk[d] )
+            throw std::logic_error( "Boundary indices out of bounds" );
+}
+
 //---------------------------------------------------------------------------//
 // Get the index space for a given combination of decomposition, entity, and
 // index types.
@@ -118,28 +154,14 @@ auto LocalGrid<MeshType>::sharedIndexSpace(
     const std::array<int, num_space_dim>& off_ijk, const int halo_width ) const
     -> IndexSpace<num_space_dim>
 {
-    // If we got the default halo width of -1 this means we want to use the
-    // default of the entire halo.
-    int hw = ( -1 == halo_width ) ? _halo_cell_width : halo_width;
-
-    // Check that the offsets are valid.
-    for ( std::size_t d = 0; d < num_space_dim; ++d )
-        if ( off_ijk[d] < -1 || 1 < off_ijk[d] )
-            throw std::logic_error( "Neighbor indices out of bounds" );
-
-    // Check that the requested halo width is valid.
-    if ( hw > _halo_cell_width )
-        throw std::logic_error(
-            "Requested halo width larger than local grid halo" );
+    auto hw = setupHaloWidth( halo_width );
+    checkOffsets( off_ijk );
 
     // Check to see if this is a valid neighbor. If not, return a shared space
     // of size 0.
     if ( neighborRank( off_ijk ) < 0 )
     {
-        std::array<long, num_space_dim> zero_size;
-        for ( std::size_t d = 0; d < num_space_dim; ++d )
-            zero_size[d] = 0;
-        return IndexSpace<num_space_dim>( zero_size, zero_size );
+        return zeroIndexSpace();
     }
 
     // Call the underlying implementation.
@@ -183,28 +205,14 @@ auto LocalGrid<MeshType>::boundaryIndexSpace(
     const std::array<int, num_space_dim>& off_ijk, const int halo_width ) const
     -> IndexSpace<num_space_dim>
 {
-    // If we got the default halo width of -1 this means we want to use the
-    // default of the entire halo.
-    int hw = ( -1 == halo_width ) ? _halo_cell_width : halo_width;
-
-    // Check that the offsets are valid.
-    for ( std::size_t d = 0; d < num_space_dim; ++d )
-        if ( off_ijk[d] < -1 || 1 < off_ijk[d] )
-            throw std::logic_error( "Boundary indices out of bounds" );
-
-    // Check that the requested halo width is valid.
-    if ( hw > _halo_cell_width )
-        throw std::logic_error(
-            "Requested halo width larger than local grid halo" );
+    auto hw = setupHaloWidth( halo_width );
+    checkOffsets( off_ijk );
 
     // Check to see if this is not a communication neighbor. If it is, return
     // a boundary space of size 0 because there is no boundary.
     if ( neighborRank( off_ijk ) >= 0 )
     {
-        std::array<long, num_space_dim> zero_size;
-        for ( std::size_t d = 0; d < num_space_dim; ++d )
-            zero_size[d] = 0;
-        return IndexSpace<num_space_dim>( zero_size, zero_size );
+        return zeroIndexSpace();
     }
 
     return boundaryIndexSpaceImpl( t1, t2, off_ijk, hw );
@@ -233,6 +241,66 @@ LocalGrid<MeshType>::boundaryIndexSpace( DecompositionTag t1, EntityType t2,
     return boundaryIndexSpace( t1, t2, off_ijk, halo_width );
 }
 
+//---------------------------------------------------------------------------//
+// Given the relative offsets of a boundary relative to this local grid's
+// indices get the set of local entity indices associated with that periodic
+// boundary in the given decomposition. Optionally provide a halo width for the
+// shared space. This halo width must be less than or equal to the halo width of
+// the local grid. The default behavior is to use the halo width of the local
+// grid. For example, if the Own decomposition is used, the interior entities
+// that would be affected by a periodic boundary operation are provided whereas
+// if the Ghost decomposition is used the halo entities on the boundary are
+// provided.
+template <class MeshType>
+template <class DecompositionTag, class EntityType>
+auto LocalGrid<MeshType>::periodicIndexSpace(
+    DecompositionTag t1, EntityType t2,
+    const std::array<int, num_space_dim>& off_ijk, const int halo_width ) const
+    -> IndexSpace<num_space_dim>
+{
+    auto hw = setupHaloWidth( halo_width );
+    checkOffsets( off_ijk );
+
+    // Check to see if this is a periodic neighbor. If not, return a boundary
+    // space of size 0 because there is no periodic boundary.
+    bool periodic_ijk = true;
+    for ( std::size_t d = 0; d < num_space_dim; ++d )
+        if ( off_ijk[d] != 0 && !_global_grid->isPeriodic( d ) )
+            periodic_ijk = false;
+
+    if ( !periodic_ijk )
+    {
+        return zeroIndexSpace();
+    }
+
+    // Call the underlying implementation. Periodic index spaces are simply the
+    // subset of shared index spaces that are on periodic boundaries.
+    return sharedIndexSpaceImpl( t1, t2, off_ijk, hw );
+}
+
+template <class MeshType>
+template <class DecompositionTag, class EntityType, std::size_t NSD>
+std::enable_if_t<3 == NSD, IndexSpace<3>>
+LocalGrid<MeshType>::periodicIndexSpace( DecompositionTag t1, EntityType t2,
+                                         const int off_i, const int off_j,
+                                         const int off_k,
+                                         const int halo_width ) const
+{
+    std::array<int, 3> off_ijk = { off_i, off_j, off_k };
+    return periodicIndexSpace( t1, t2, off_ijk, halo_width );
+}
+
+template <class MeshType>
+template <class DecompositionTag, class EntityType, std::size_t NSD>
+std::enable_if_t<2 == NSD, IndexSpace<2>>
+LocalGrid<MeshType>::periodicIndexSpace( DecompositionTag t1, EntityType t2,
+                                         const int off_i, const int off_j,
+                                         const int halo_width ) const
+{
+    std::array<int, 2> off_ijk = { off_i, off_j };
+    return periodicIndexSpace( t1, t2, off_ijk, halo_width );
+}
+
 //---------------------------------------------------------------------------//
 // Get the local index space of the owned cells.
 template <class MeshType>