diff --git a/example/grid_tutorial/06_local_grid/local_grid_example.cpp b/example/grid_tutorial/06_local_grid/local_grid_example.cpp index 20af3a0db..cee5352ed 100644 --- a/example/grid_tutorial/06_local_grid/local_grid_example.cpp +++ b/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. diff --git a/grid/src/Cabana_Grid_LocalGrid.hpp b/grid/src/Cabana_Grid_LocalGrid.hpp index e086dfdae..39347d9fe 100644 --- a/grid/src/Cabana_Grid_LocalGrid.hpp +++ b/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 + IndexSpace + periodicIndexSpace( DecompositionTag t1, EntityType t2, + const std::array& 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 + 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 + 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& off_ijk ) const; + template auto getBound( OwnedIndexSpace owned_space, const int upper_lower, const std::array& off_ijk, diff --git a/grid/src/Cabana_Grid_LocalGrid_impl.hpp b/grid/src/Cabana_Grid_LocalGrid_impl.hpp index cad4031d5..d8b54017e 100644 --- a/grid/src/Cabana_Grid_LocalGrid_impl.hpp +++ b/grid/src/Cabana_Grid_LocalGrid_impl.hpp @@ -92,6 +92,42 @@ LocalGrid::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 +auto LocalGrid::zeroIndexSpace() const +{ + std::array zero_size; + for ( std::size_t d = 0; d < num_space_dim; ++d ) + zero_size[d] = 0; + return IndexSpace( zero_size, zero_size ); +} + +template +auto LocalGrid::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 +void LocalGrid::checkOffsets( + const std::array& 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::sharedIndexSpace( const std::array& off_ijk, const int halo_width ) const -> IndexSpace { - // 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 zero_size; - for ( std::size_t d = 0; d < num_space_dim; ++d ) - zero_size[d] = 0; - return IndexSpace( zero_size, zero_size ); + return zeroIndexSpace(); } // Call the underlying implementation. @@ -183,28 +205,14 @@ auto LocalGrid::boundaryIndexSpace( const std::array& off_ijk, const int halo_width ) const -> IndexSpace { - // 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 zero_size; - for ( std::size_t d = 0; d < num_space_dim; ++d ) - zero_size[d] = 0; - return IndexSpace( zero_size, zero_size ); + return zeroIndexSpace(); } return boundaryIndexSpaceImpl( t1, t2, off_ijk, hw ); @@ -233,6 +241,66 @@ LocalGrid::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 +template +auto LocalGrid::periodicIndexSpace( + DecompositionTag t1, EntityType t2, + const std::array& off_ijk, const int halo_width ) const + -> IndexSpace +{ + 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 +template +std::enable_if_t<3 == NSD, IndexSpace<3>> +LocalGrid::periodicIndexSpace( DecompositionTag t1, EntityType t2, + const int off_i, const int off_j, + const int off_k, + const int halo_width ) const +{ + std::array off_ijk = { off_i, off_j, off_k }; + return periodicIndexSpace( t1, t2, off_ijk, halo_width ); +} + +template +template +std::enable_if_t<2 == NSD, IndexSpace<2>> +LocalGrid::periodicIndexSpace( DecompositionTag t1, EntityType t2, + const int off_i, const int off_j, + const int halo_width ) const +{ + std::array off_ijk = { off_i, off_j }; + return periodicIndexSpace( t1, t2, off_ijk, halo_width ); +} + //---------------------------------------------------------------------------// // Get the local index space of the owned cells. template