From d657b0fbee5e92c2a5418930ea520e88407a0312 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 9 Aug 2024 09:38:30 -0400 Subject: [PATCH 1/3] fixup: local grid index space helper functions --- grid/src/Cabana_Grid_LocalGrid.hpp | 4 ++ grid/src/Cabana_Grid_LocalGrid_impl.hpp | 73 +++++++++++++------------ 2 files changed, 43 insertions(+), 34 deletions(-) diff --git a/grid/src/Cabana_Grid_LocalGrid.hpp b/grid/src/Cabana_Grid_LocalGrid.hpp index e086dfdae..bc118295c 100644 --- a/grid/src/Cabana_Grid_LocalGrid.hpp +++ b/grid/src/Cabana_Grid_LocalGrid.hpp @@ -253,6 +253,10 @@ class LocalGrid private: // Helper functions + zeroIndexSpace(); + setupHaloWidthImpl( const int halo_width ); + checkOffsets( const std::array& off_ijk ); + 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..584fd1097 100644 --- a/grid/src/Cabana_Grid_LocalGrid_impl.hpp +++ b/grid/src/Cabana_Grid_LocalGrid_impl.hpp @@ -92,6 +92,39 @@ 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. +auto LocalGrid::zeroIndexSpace() +{ + 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 ); +} + +auto LocalGrid::setupHaloWidth( const int halo_width ) +{ + // 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; +} + +auto LocalGrid::checkOffsets( + const std::array& off_ijk ) +{ + // 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 +151,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 +202,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 ); From c64f42de8bc932dadbae49fb673993d3cb245ba8 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 9 Aug 2024 09:40:46 -0400 Subject: [PATCH 2/3] Add periodicIndexSpace analogous to boundaryIndexSpace; subset of sharedIndexSpace on periodic boundaries --- grid/src/Cabana_Grid_LocalGrid.hpp | 76 ++++++++++++++++++++++++- grid/src/Cabana_Grid_LocalGrid_impl.hpp | 71 +++++++++++++++++++++-- 2 files changed, 140 insertions(+), 7 deletions(-) diff --git a/grid/src/Cabana_Grid_LocalGrid.hpp b/grid/src/Cabana_Grid_LocalGrid.hpp index bc118295c..39347d9fe 100644 --- a/grid/src/Cabana_Grid_LocalGrid.hpp +++ b/grid/src/Cabana_Grid_LocalGrid.hpp @@ -251,11 +251,81 @@ 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 - zeroIndexSpace(); - setupHaloWidthImpl( const int halo_width ); - checkOffsets( const std::array& off_ijk ); + 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, diff --git a/grid/src/Cabana_Grid_LocalGrid_impl.hpp b/grid/src/Cabana_Grid_LocalGrid_impl.hpp index 584fd1097..d8b54017e 100644 --- a/grid/src/Cabana_Grid_LocalGrid_impl.hpp +++ b/grid/src/Cabana_Grid_LocalGrid_impl.hpp @@ -95,7 +95,8 @@ LocalGrid::neighborRank( const int off_i, const int off_j ) const //---------------------------------------------------------------------------// // Get the index space for a given combination of decomposition, entity, and // index types. -auto LocalGrid::zeroIndexSpace() +template +auto LocalGrid::zeroIndexSpace() const { std::array zero_size; for ( std::size_t d = 0; d < num_space_dim; ++d ) @@ -103,7 +104,8 @@ auto LocalGrid::zeroIndexSpace() return IndexSpace( zero_size, zero_size ); } -auto LocalGrid::setupHaloWidth( const int halo_width ) +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. @@ -116,8 +118,9 @@ auto LocalGrid::setupHaloWidth( const int halo_width ) return hw; } -auto LocalGrid::checkOffsets( - const std::array& off_ijk ) +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 ) @@ -238,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 From fab6489238bc5daaa4e4f5ec52d8dd60e24b73f5 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Thu, 19 Sep 2024 10:26:36 -0400 Subject: [PATCH 3/3] Add periodic index space to local grid example --- .../06_local_grid/local_grid_example.cpp | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) 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.