Cabana_Grid_SparseIndexSpace.hpp

Go to the documentation of this file.

/****************************************************************************
 * Copyright (c) 2018-2023 by the Cabana authors                            *
 * All rights reserved.                                                     *
 *                                                                          *
 * This file is part of the Cabana library. Cabana is distributed under a   *
 * BSD 3-clause license. For the licensing terms see the LICENSE file in    *
 * the top-level directory.                                                 *
 *                                                                          *
 * SPDX-License-Identifier: BSD-3-Clause                                    *
 ****************************************************************************/

#ifndef CABANA_GRID_SPARSE_INDEXSPACE_HPP
#define CABANA_GRID_SPARSE_INDEXSPACE_HPP
 
#ifdef KOKKOS_ENABLE_SYCL
#error Not supported with the Kokkos SYCL backend (device recursion)
#endif
 
#include <Cabana_Grid_GlobalMesh.hpp>
#include <Cabana_Grid_IndexSpace.hpp>
 
#include <Kokkos_Core.hpp>
#include <Kokkos_UnorderedMap.hpp>
 
#include <array>
#include <memory>
#include <string>
 
//---------------------------------------------------------------------------//
// Naming convention:
//  Tile ID / Tile ijk = 3D tile indexing, i.e., (i, j, k)
//  Tile No. = 1D number indicating the tile position in the allocated array
//  Tile key / Hash key = 1D number, computed with the given indexing method and
//  used as hash key
//---------------------------------------------------------------------------//
 
namespace Cabana
{
namespace Grid
{
 
//---------------------------------------------------------------------------//
enum class HashTypes : unsigned char
{
    Naive = 0, // Lexicographical Order
    Morton = 1 // Morton Curve
};
 
//---------------------------------------------------------------------------//
// bit operations for Morton Code computing
//---------------------------------------------------------------------------//
 
//---------------------------------------------------------------------------//
template <typename Integer>
KOKKOS_INLINE_FUNCTION constexpr Integer bitLength( Integer input_int ) noexcept
{
    if ( input_int )
        return bitLength( input_int >> 1 ) + static_cast<Integer>( 1 );
    else
        return 0;
}
 
//---------------------------------------------------------------------------//
template <typename Integer>
KOKKOS_INLINE_FUNCTION constexpr Integer bitCount( Integer input_int ) noexcept
{
    return bitLength( input_int - 1 );
}
 
//---------------------------------------------------------------------------//
template <typename Integer>
KOKKOS_INLINE_FUNCTION constexpr Integer
binaryReverse( Integer input_int, char loc = sizeof( Integer ) * 8 - 1 )
{
    if ( input_int == 0 )
        return 0;
    return ( ( input_int & 1 ) << loc ) |
           binaryReverse( input_int >> 1, loc - 1 );
}
 
//---------------------------------------------------------------------------//
template <typename Integer>
KOKKOS_INLINE_FUNCTION constexpr unsigned countLeadingZeros( Integer input_int )
{
    unsigned res{ 0 };
    input_int = binaryReverse( input_int );
    if ( input_int == 0 )
        return sizeof( Integer ) * 8;
    while ( ( input_int & 1 ) == 0 )
        res++, input_int >>= 1;
    return res;
}
 
//---------------------------------------------------------------------------//
KOKKOS_INLINE_FUNCTION
constexpr int bitPack( const uint64_t mask, const uint64_t data )
{
    uint64_t slresult = 0;
    uint64_t& ulresult{ slresult };
    uint64_t uldata = data;
    int count = 0;
    ulresult = 0;
 
    uint64_t rmask = binaryReverse( mask );
    unsigned char lz{ 0 };
 
    while ( rmask )
    {
        lz = countLeadingZeros( rmask );
        uldata >>= lz;
        ulresult <<= 1;
        count++;
        ulresult |= ( uldata & 1 );
        uldata >>= 1;
        rmask <<= lz + 1;
    }
    ulresult <<= 64 - count; // 64 bits (maybe not use a constant 64 ...?)
    ulresult = binaryReverse( ulresult );
    return (int)slresult;
}
 
//---------------------------------------------------------------------------//
KOKKOS_INLINE_FUNCTION
constexpr uint64_t bitSpread( const uint64_t mask, const int data )
{
    uint64_t rmask = binaryReverse( mask );
    int dat = data;
    uint64_t result = 0;
    unsigned char lz{ 0 };
    while ( rmask )
    {
        lz = countLeadingZeros( rmask ) + 1;
        result = result << lz | ( dat & 1 );
        dat >>= 1, rmask <<= lz;
    }
    result = binaryReverse( result ) >> countLeadingZeros( mask );
    return result;
}
 
//---------------------------------------------------------------------------//
// Tile Hash Key (hash) <=> TileID.
//---------------------------------------------------------------------------//
template <typename Key, HashTypes HashT>
struct TileID2HashKey;
 
template <typename Key, HashTypes HashT>
struct HashKey2TileID;

template <typename Key>
struct TileID2HashKey<Key, HashTypes::Naive>
{
    using tid_to_key_type = TileID2HashKey<Key, HashTypes::Naive>;
    static constexpr HashTypes hash_type = HashTypes::Naive;

    TileID2HashKey( std::initializer_list<int> tnum )
    {
        std::copy( tnum.begin(), tnum.end(), _tile_num.data() );
        _tile_num_j_mul_k = _tile_num[1] * _tile_num[2];
    }
    TileID2HashKey( int i, int j, int k )
    {
        _tile_num[0] = i;
        _tile_num[1] = j;
        _tile_num[2] = k;
        _tile_num_j_mul_k = j * k;
    }

    KOKKOS_INLINE_FUNCTION
    Key operator()( int tile_i, int tile_j, int tile_k ) const
    {
        return tile_i * _tile_num_j_mul_k + tile_j * _tile_num[2] + tile_k;
    }
 
  private:
    // The index bounds of the tiles on the current MPI rank
    Kokkos::Array<int, 3> _tile_num;
    int _tile_num_j_mul_k;
};

template <typename Key>
struct TileID2HashKey<Key, HashTypes::Morton>
{
    using tid_to_key_type = TileID2HashKey<Key, HashTypes::Morton>;
    static constexpr HashTypes hash_type = HashTypes::Morton;

    TileID2HashKey( std::initializer_list<int> tnum )
    {
        std::copy( tnum.begin(), tnum.end(), _tile_num.data() );
    }
    TileID2HashKey( int i, int j, int k )
    {
        _tile_num[0] = i;
        _tile_num[1] = j;
        _tile_num[2] = k;
    }
    enum : uint64_t
    { // hand-coded now, can be improved by iterative func
        page_kmask = ( 0x9249249249249249UL ),
        page_jmask = ( 0x2492492492492492UL ),
        page_imask = ( 0x4924924924924924UL )
    };

    KOKKOS_INLINE_FUNCTION
    Key operator()( int tile_i, int tile_j, int tile_k ) const
    {
        return bitSpread( page_kmask, tile_k ) |
               bitSpread( page_jmask, tile_j ) |
               bitSpread( page_imask, tile_i );
    }
 
  private:
    // The index bounds of the tiles on the current MPI rank
    Kokkos::Array<int, 3> _tile_num;
};

template <typename Key>
struct HashKey2TileID<Key, HashTypes::Naive>
{
    using key_to_tid_type = HashKey2TileID<Key, HashTypes::Naive>;
    static constexpr HashTypes hash_type = HashTypes::Naive;

    HashKey2TileID( std::initializer_list<int> tnum )
        : _tile_num()
    {
        std::copy( tnum.begin(), tnum.end(), _tile_num.data() );
    }
    HashKey2TileID( int i, int j, int k )
    {
        _tile_num[0] = i;
        _tile_num[1] = j;
        _tile_num[2] = k;
    }

    KOKKOS_INLINE_FUNCTION
    void operator()( Key tile_key, int& tile_i, int& tile_j, int& tile_k ) const
    {
        tile_k = tile_key % _tile_num[2];
        tile_j = static_cast<Key>( tile_key / _tile_num[2] ) % _tile_num[1];
        tile_i = static_cast<Key>( tile_key / _tile_num[2] / _tile_num[1] ) %
                 _tile_num[0];
    }
 
  private:
    // The index bounds of the tiles on the current MPI rank
    Kokkos::Array<int, 3> _tile_num;
};

template <typename Key>
struct HashKey2TileID<Key, HashTypes::Morton>
{
    using key_to_tid_type = HashKey2TileID<Key, HashTypes::Morton>;
    static constexpr HashTypes hash_type = HashTypes::Morton;

    HashKey2TileID( std::initializer_list<int> tnum )
    {
        std::copy( tnum.begin(), tnum.end(), _tile_num.data() );
    }
    HashKey2TileID( int i, int j, int k )
    {
        _tile_num[0] = i;
        _tile_num[1] = j;
        _tile_num[2] = k;
    }
    enum : uint64_t
    { // hand-coded now, can be improved by iterative func
        page_kmask = ( 0x9249249249249249UL ),
        page_jmask = ( 0x2492492492492492UL ),
        page_imask = ( 0x4924924924924924UL )
    };

    KOKKOS_INLINE_FUNCTION
    void operator()( Key tile_key, int& tile_i, int& tile_j, int& tile_k ) const
    {
        tile_k = bitPack( page_kmask, tile_key );
        tile_j = bitPack( page_jmask, tile_key );
        tile_i = bitPack( page_imask, tile_key );
    }
 
  private:
    // The index bounds of the tiles on the current MPI rank
    Kokkos::Array<int, 3> _tile_num;
};
 
//---------------------------------------------------------------------------//
// Hierarchical index spaces
// SparseMap <- BlockMap <- TileMap
// Naming:
//      Block = MPI Rank
//      Tile = Sub-block with several cells
//      Cell = Basic grid unit
//---------------------------------------------------------------------------//

template <typename MemorySpace, unsigned long long CBits,
          unsigned long long CNumPerDim, unsigned long long CNumPerTile,
          HashTypes Hash, typename Key, typename Value>
class BlockMap;

template <int CBits, int CNumPerDim, int CNumPerTile>
class TileMap;

template <typename MemorySpace, unsigned long long CellPerTileDim = 4,
          HashTypes Hash = HashTypes::Naive, typename Key = uint64_t,
          typename Value = uint64_t>
class SparseMap
{
  public:
    static constexpr int rank = 3;
    static constexpr unsigned long long cell_bits_per_tile_dim =
        bitCount( CellPerTileDim );
    static constexpr unsigned long long cell_num_per_tile_dim =
        1 << cell_bits_per_tile_dim;
    static constexpr unsigned long long cell_mask_per_tile_dim =
        ( 1 << cell_bits_per_tile_dim ) - 1;
    static constexpr unsigned long long cell_bits_per_tile =
        cell_bits_per_tile_dim + cell_bits_per_tile_dim +
        cell_bits_per_tile_dim;
    static constexpr unsigned long long cell_mask_per_tile =
        ( 1 << cell_bits_per_tile ) - 1;
    static constexpr unsigned long long cell_num_per_tile =
        cell_num_per_tile_dim * cell_num_per_tile_dim * cell_num_per_tile_dim;
    using key_type = Key;
    using value_type = Value;
    static constexpr HashTypes hash_type = Hash;
    using memory_space = MemorySpace;

    SparseMap( const std::array<int, rank> size,
               const unsigned int pre_alloc_size )
        : _block_id_space( size[0] >> cell_bits_per_tile_dim,
                           size[1] >> cell_bits_per_tile_dim,
                           size[2] >> cell_bits_per_tile_dim,
                           1 << bitCount( pre_alloc_size ) )
    {
        std::fill( _min.data(), _min.data() + rank, 0 );
        std::copy( size.begin(), size.end(), _max.data() );
    }

    KOKKOS_INLINE_FUNCTION
    void insertCell( int cell_i, int cell_j, int cell_k ) const
    {
        insertTile( cell_i >> cell_bits_per_tile_dim,
                    cell_j >> cell_bits_per_tile_dim,
                    cell_k >> cell_bits_per_tile_dim );
    }

    KOKKOS_INLINE_FUNCTION
    void insertTile( int tile_i, int tile_j, int tile_k ) const
    {
        _block_id_space.insert( tile_i, tile_j, tile_k );
    }

    KOKKOS_INLINE_FUNCTION
    value_type queryTile( int cell_i, int cell_j, int cell_k ) const
    {
        // query the tile ID
        auto tile_id = _block_id_space.find( cell_i >> cell_bits_per_tile_dim,
                                             cell_j >> cell_bits_per_tile_dim,
                                             cell_k >> cell_bits_per_tile_dim );
        return tile_id;
    }

    KOKKOS_INLINE_FUNCTION
    value_type queryTileFromTileId( const int tile_i, const int tile_j,
                                    const int tile_k ) const
    {
        // query the tile ID
        return _block_id_space.find( tile_i, tile_j, tile_k );
    }

    KOKKOS_INLINE_FUNCTION
    value_type queryTileFromTileKey( const key_type tile_key ) const
    {
        return _block_id_space.findKey( tile_key );
    }

    KOKKOS_INLINE_FUNCTION
    value_type queryCell( int cell_i, int cell_j, int cell_k ) const
    {
        // query the tile ID.
        auto tile_id = queryTile( cell_i, cell_j, cell_k );
        auto cell_id = cell_local_id( cell_i, cell_j, cell_k );
        return static_cast<value_type>( ( tile_id << cell_bits_per_tile ) |
                                        cell_id );
    }

    void clear() { _block_id_space.clear(); }

    bool reserve( const value_type capacity )
    {
        return _block_id_space.reserve( capacity );
    }

    KOKKOS_INLINE_FUNCTION
    uint32_t capacity() const { return _block_id_space.capacity(); }

    value_type size() const { return _block_id_space.validTileNumHost(); }

    inline value_type sizeTile() const { return size(); }

    inline value_type sizeCell() const
    {
        return sizeTile() << cell_bits_per_tile;
    }

    template <typename T>
    inline value_type reservedTileSize( T factor ) const
    {
        return _block_id_space.reservedTileNumHost( factor );
    }

    template <typename T>
    inline value_type reservedCellSize( T factor ) const
    {
        return reservedTileSize( factor ) << cell_bits_per_tile;
    }

    KOKKOS_INLINE_FUNCTION
    bool valid_at( uint32_t index ) const
    {
        return _block_id_space.valid_at( index );
    }

    KOKKOS_INLINE_FUNCTION
    bool isValidKey( key_type key ) const
    {
        return _block_id_space.isValidKey( key );
    }

    KOKKOS_INLINE_FUNCTION
    key_type key_at( uint32_t index ) const
    {
        return _block_id_space.key_at( index );
    }

    KOKKOS_INLINE_FUNCTION
    value_type value_at( uint32_t index ) const
    {
        return _block_id_space.value_at( index );
    }

    KOKKOS_INLINE_FUNCTION
    void key2ijk( key_type& key, int& tile_i, int& tile_j, int& tile_k ) const
    {
        return _block_id_space.key2ijk( key, tile_i, tile_j, tile_k );
    }

    KOKKOS_INLINE_FUNCTION
    key_type ijk2key( int tile_i, int tile_j, int tile_k ) const
    {
        return _block_id_space.ijk2key( tile_i, tile_j, tile_k );
    }

    KOKKOS_INLINE_FUNCTION
    value_type cell_local_id( const int cell_i, const int cell_j,
                              const int cell_k ) const
    {
        return _tile_id_space.coordToOffset( cell_i & cell_mask_per_tile_dim,
                                             cell_j & cell_mask_per_tile_dim,
                                             cell_k & cell_mask_per_tile_dim );
    }
 
  private:
    BlockMap<MemorySpace, cell_bits_per_tile_dim, cell_num_per_tile_dim,
             cell_num_per_tile, hash_type, key_type, value_type>
        _block_id_space;
    TileMap<cell_bits_per_tile_dim, cell_num_per_tile_dim, cell_num_per_tile>
        _tile_id_space;
    Kokkos::Array<int, rank> _min;
    Kokkos::Array<int, rank> _max;
};
 
//---------------------------------------------------------------------------//
template <typename MemorySpace, class Scalar,
          unsigned long long CellPerTileDim = 4,
          HashTypes Hash = HashTypes::Naive, typename Key = uint64_t,
          typename Value = uint64_t>
SparseMap<MemorySpace, CellPerTileDim, Hash, Key, Value> createSparseMap(
    const std::shared_ptr<GlobalMesh<SparseMesh<Scalar>>>& global_mesh,
    int pre_alloc_size )
{
    return SparseMap<MemorySpace, CellPerTileDim, Hash, Key, Value>(
        { global_mesh->globalNumCell( Dim::I ),
          global_mesh->globalNumCell( Dim::J ),
          global_mesh->globalNumCell( Dim::K ) },
        pre_alloc_size );
}
 
//---------------------------------------------------------------------------//
template <typename MemorySpace, unsigned long long CBits,
          unsigned long long CNumPerDim, unsigned long long CNumPerTile,
          HashTypes Hash, typename Key, typename Value>
class BlockMap
{
  public:
    static constexpr unsigned long long cell_bits_per_tile_dim = CBits;
    static constexpr unsigned long long cell_num_per_tile_dim = CNumPerDim;
    static constexpr unsigned long long cell_num_per_tile = CNumPerTile;
    using key_type = Key;
    using value_type = Value;
    static constexpr HashTypes hash_type = Hash;
    using bis_Type =
        BlockMap<MemorySpace, cell_bits_per_tile_dim, cell_num_per_tile_dim,
                 cell_num_per_tile, hash_type, key_type, value_type>; // itself

    BlockMap( const int size_x, const int size_y, const int size_z,
              const value_type pre_alloc_size )
        : _tile_table_info( "hash_table_info" )
        , _tile_table( size_x * size_y * size_z )
        , _op_ijk2key( size_x, size_y, size_z )
        , _op_key2ijk( size_x, size_y, size_z )
    {
        // hash table related init
        auto tile_table_info_mirror =
            Kokkos::create_mirror_view( Kokkos::HostSpace(), _tile_table_info );
        tile_table_info_mirror( 0 ) = 0;
        tile_table_info_mirror( 1 ) = pre_alloc_size;
        Kokkos::deep_copy( _tile_table_info, tile_table_info_mirror );
 
        // size related init
        _block_size[0] = size_x;
        _block_size[1] = size_y;
        _block_size[2] = size_z;
    }

    void clear()
    {
        _tile_table.clear(); // clear hash table
        auto tile_table_info_mirror = Kokkos::create_mirror_view_and_copy(
            Kokkos::HostSpace(), _tile_table_info );
        tile_table_info_mirror( 0 ) = 0;
        Kokkos::deep_copy( _tile_table_info,
                           tile_table_info_mirror ); // clear valid size record
    }

    bool reserve( const value_type capacity )
    {
        return _tile_table.rehash( capacity );
    }

    KOKKOS_INLINE_FUNCTION
    uint32_t capacity() const
    {
        return _tile_table.capacity(); // hash_table capacity
    }

    KOKKOS_INLINE_FUNCTION
    bool valid_at( uint32_t index ) const
    {
        return _tile_table.valid_at( index );
    }

    KOKKOS_INLINE_FUNCTION
    bool isValidKey( key_type key ) const
    {
        return !( _tile_table.find( key ) == ~static_cast<uint32_t>( 0 ) );
    }

    KOKKOS_INLINE_FUNCTION
    key_type key_at( uint32_t index ) const
    {
        return _tile_table.key_at( index );
    }

    KOKKOS_INLINE_FUNCTION
    key_type value_at( uint32_t index ) const
    {
        return _tile_table.value_at( index );
    }

    KOKKOS_INLINE_FUNCTION
    value_type validTileNumDev() const { return _tile_table_info( 0 ); }

    value_type validTileNumHost() const
    {
        auto tile_table_info_mirror = Kokkos::create_mirror_view_and_copy(
            Kokkos::HostSpace(), _tile_table_info );
        return tile_table_info_mirror( 0 );
    }

    template <typename T>
    inline value_type reservedTileNumHost( T factor ) const
    {
        auto tile_table_info_mirror = Kokkos::create_mirror_view_and_copy(
            Kokkos::HostSpace(), _tile_table_info );
        value_type num = tile_table_info_mirror( 0 ) * factor;
        tile_table_info_mirror( 1 ) = num;
        Kokkos::deep_copy( _tile_table_info, tile_table_info_mirror );
        return num;
    }

    KOKKOS_INLINE_FUNCTION
    void insert( int tile_i, int tile_j, int tile_k ) const
    {
        // Note: reallocation will be performed after insert
        // tile ijk => tile key
        key_type tile_key = ijk2key( tile_i, tile_j, tile_k );
        // Try to insert the tile key into the map. Give it a dummy value for
        // now.
        auto insert_result = _tile_table.insert( tile_key, 0 );
        // If the tile key was actually inserted, atomically increment the
        // counter Only threads that actually did a successful insert will call
        // it
        if ( !insert_result.existing() )
            _tile_table.value_at( insert_result.index() ) =
                Kokkos::atomic_fetch_add( &( _tile_table_info( 0 ) ), 1 );
    }

    KOKKOS_INLINE_FUNCTION
    value_type find( int tile_i, int tile_j, int tile_k ) const
    {
        return _tile_table.value_at(
            _tile_table.find( ijk2key( tile_i, tile_j, tile_k ) ) );
    }

    KOKKOS_INLINE_FUNCTION
    value_type findKey( key_type tile_key ) const
    {
        return _tile_table.value_at( _tile_table.find( tile_key ) );
    }

    KOKKOS_INLINE_FUNCTION
    key_type ijk2key( int tile_i, int tile_j, int tile_k ) const
    {
        return _op_ijk2key( tile_i, tile_j, tile_k );
    }

    KOKKOS_INLINE_FUNCTION
    void key2ijk( key_type& key, int& tile_i, int& tile_j, int& tile_k ) const
    {
        return _op_key2ijk( key, tile_i, tile_j, tile_k );
    }
 
  private:
    Kokkos::View<int[2], MemorySpace> _tile_table_info;
    Kokkos::Array<int, 3> _block_size;
    Kokkos::UnorderedMap<key_type, value_type, MemorySpace> _tile_table;
    TileID2HashKey<key_type, hash_type> _op_ijk2key;
    HashKey2TileID<key_type, hash_type> _op_key2ijk;
};
 
//---------------------------------------------------------------------------//
template <int CBits, int CNumPerDim, int CNumPerTile>
class TileMap
{
  public:
    static constexpr int cell_bits_per_tile_dim = CBits;
    static constexpr int cell_num_per_tile_dim = CNumPerDim;
    static constexpr int cell_num_per_tile = CNumPerTile;
    static constexpr int rank = 3;


    template <typename... Coords>
    KOKKOS_INLINE_FUNCTION static constexpr auto
    coordToOffset( Coords&&... coords ) -> uint64_t
    {
        return fromCoord<Coord2OffsetDim>( std::forward<Coords>( coords )... );
    }

    template <typename Key, typename... Coords>
    KOKKOS_INLINE_FUNCTION static constexpr void
    offsetToCoord( Key&& key, Coords&... coords )
    {
        toCoord<Offset2CoordDim>( std::forward<Key>( key ), coords... );
    }
 
  private:

    struct Coord2OffsetDim
    {
        template <typename Coord>
        KOKKOS_INLINE_FUNCTION constexpr auto operator()( int dim_no,
                                                          Coord&& i )
            -> uint64_t
        {
            uint64_t result = static_cast<uint64_t>( i );
            for ( int i = 0; i < dim_no; i++ )
                result *= cell_num_per_tile_dim;
            return result;
        }
    };

    struct Offset2CoordDim
    {
        template <typename Coord, typename Key>
        KOKKOS_INLINE_FUNCTION constexpr auto operator()( Coord& i,
                                                          Key&& offset )
            -> uint64_t
        {
            i = offset % cell_num_per_tile_dim;
            return ( offset / cell_num_per_tile_dim );
        }
    };

    template <typename Func, typename... Coords>
    KOKKOS_INLINE_FUNCTION static constexpr auto fromCoord( Coords&&... coords )
        -> uint64_t
    {
        static_assert( sizeof...( Coords ) == rank,
                       "Dimension of coordinate mismatch" );
        using integer = std::common_type_t<Coords...>;
        static_assert( std::is_integral<integer>::value,
                       "Coordinate should be integral type" );
        return fromCoordImpl<Func>( 0, coords... );
    }

    template <typename Func, typename Key, typename... Coords>
    KOKKOS_INLINE_FUNCTION static constexpr void toCoord( Key&& key,
                                                          Coords&... coords )
    {
        static_assert( sizeof...( Coords ) == rank,
                       "Dimension of coordinate mismatch" );
        using integer = std::common_type_t<Coords...>;
        static_assert( std::is_integral<integer>::value,
                       "Coordinate should be integral type" );
        return toCoordImpl<Func>( 0, std::forward<Key>( key ), coords... );
    }

    template <typename Func, typename Coord>
    KOKKOS_INLINE_FUNCTION static constexpr auto fromCoordImpl( int&& dim_no,
                                                                Coord&& i )
    {
        return Func()( dim_no, std::forward<Coord>( i ) );
    }

    template <typename Func, typename Coord, typename... Coords>
    KOKKOS_INLINE_FUNCTION static constexpr auto
    fromCoordImpl( int&& dim_no, Coord&& i, Coords&&... is )
    {
        auto result = Func()( dim_no, std::forward<Coord>( i ) );
        if ( dim_no + 1 < rank )
            result += fromCoordImpl<Func>( dim_no + 1,
                                           std::forward<Coords>( is )... );
        return result;
    }

    template <typename Func, typename Key, typename Coord>
    KOKKOS_INLINE_FUNCTION static constexpr void
    toCoordImpl( int&& dim_no, Key&& key, Coord& i ) noexcept
    {
        if ( dim_no < rank )
            Func()( i, std::forward<Key>( key ) );
    }

    template <typename Func, typename Key, typename Coord, typename... Coords>
    KOKKOS_INLINE_FUNCTION static constexpr void
    toCoordImpl( int&& dim_no, Key&& key, Coord& i, Coords&... is ) noexcept
    {
        auto new_key = Func()( i, std::forward<Key>( key ) );
        if ( dim_no + 1 < rank )
            toCoordImpl<Func>( dim_no + 1, new_key, is... );
    }
};
 
template <long N>
class IndexSpace;
 
//---------------------------------------------------------------------------//
template <std::size_t N, unsigned long long cellBitsPerTileDim,
          typename std::enable_if_t<( N == 3 ), bool> = true>
class TileIndexSpace : public IndexSpace<N>
{
  public:
    static constexpr std::size_t Rank = N;
    static constexpr unsigned long long cell_bits_per_tile_dim =
        cellBitsPerTileDim;
    static constexpr unsigned long long cell_num_per_tile =
        1 << ( cell_bits_per_tile_dim * Rank );

    TileIndexSpace()
        : IndexSpace<N>()
    {
    }

    TileIndexSpace( const std::initializer_list<long>& size )
        : IndexSpace<N>( size )
    {
    }

    TileIndexSpace( const std::initializer_list<long>& min,
                    const std::initializer_list<long>& max )
        : IndexSpace<N>( min, max )
    {
    }

    template <typename... Params>
    TileIndexSpace( Params&&... pars )
        : IndexSpace<N>( std::forward<Params>( pars )... )
    {
    }

    KOKKOS_INLINE_FUNCTION
    long min( const long dim ) const { return IndexSpace<N>::min( dim ); }

    KOKKOS_INLINE_FUNCTION
    long max( const long dim ) const { return IndexSpace<N>::max( dim ); }

    KOKKOS_INLINE_FUNCTION
    Kokkos::Array<long, Rank> min() const { return IndexSpace<N>::min(); }

    KOKKOS_INLINE_FUNCTION
    Kokkos::Array<long, Rank> max() const { return IndexSpace<N>::max(); }

    KOKKOS_FORCEINLINE_FUNCTION
    long sizeTile() const { return IndexSpace<N>::size(); }

    KOKKOS_FORCEINLINE_FUNCTION
    long sizeCell() const { return ( sizeTile() * cell_num_per_tile ); }

    template <int NSD = N>
    KOKKOS_FORCEINLINE_FUNCTION std::enable_if_t<3 == NSD, bool>
    tileInRange( const long tile_i, const long tile_j, const long tile_k ) const
    {
        long index[N] = { tile_i, tile_j, tile_k };
        return IndexSpace<N>::inRange( index );
    }

    template <int NSD = N>
    KOKKOS_FORCEINLINE_FUNCTION std::enable_if_t<3 == NSD, bool>
    cellInRange( const long cell_i, const long cell_j, const long cell_k ) const
    {
        return tileInRange( cell_i >> cell_bits_per_tile_dim,
                            cell_j >> cell_bits_per_tile_dim,
                            cell_k >> cell_bits_per_tile_dim );
    }
};
 
} // namespace Grid
} // namespace Cabana
 
#endif