Cabana_Grid_Array.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_ARRAY_HPP
#define CABANA_GRID_ARRAY_HPP
 
#include <Cabana_Grid_IndexSpace.hpp>
#include <Cabana_Grid_LocalGrid.hpp>
#include <Cabana_Grid_MpiTraits.hpp>
#include <Cabana_Grid_Types.hpp>
 
#include <Kokkos_Core.hpp>
 
#include <cmath>
#include <memory>
#include <type_traits>
#include <vector>
 
#include <mpi.h>
 
namespace Cabana
{
namespace Grid
{
//---------------------------------------------------------------------------//
template <class EntityType, class MeshType>
class ArrayLayout
{
  public:
    using entity_type = EntityType;

    using mesh_type = MeshType;

    static constexpr std::size_t num_space_dim = mesh_type::num_space_dim;

    ArrayLayout( const std::shared_ptr<LocalGrid<MeshType>>& local_grid,
                 const int dofs_per_entity )
        : _local_grid( local_grid )
        , _dofs_per_entity( dofs_per_entity )
    {
    }

    const std::shared_ptr<LocalGrid<MeshType>> localGrid() const
    {
        return _local_grid;
    }

    int dofsPerEntity() const { return _dofs_per_entity; }

    template <class DecompositionTag, class IndexType>
    IndexSpace<num_space_dim + 1>
    indexSpace( DecompositionTag decomposition_tag, IndexType index_type ) const
    {
        return appendDimension( _local_grid->indexSpace( decomposition_tag,
                                                         EntityType(),
                                                         index_type ),
                                _dofs_per_entity );
    }

    template <class DecompositionTag>
    IndexSpace<num_space_dim + 1>
    sharedIndexSpace( DecompositionTag decomposition_tag,
                      const std::array<int, num_space_dim>& off_ijk,
                      const int halo_width = -1 ) const
    {
        return appendDimension(
            _local_grid->sharedIndexSpace( decomposition_tag, EntityType(),
                                           off_ijk, halo_width ),
            _dofs_per_entity );
    }

    template <class DecompositionTag, std::size_t NSD = num_space_dim>
    std::enable_if_t<3 == NSD, IndexSpace<4>>
    sharedIndexSpace( DecompositionTag decomposition_tag, const int off_i,
                      const int off_j, const int off_k,
                      const int halo_width = -1 ) const
    {
        std::array<int, 3> off_ijk = { off_i, off_j, off_k };
        return sharedIndexSpace( decomposition_tag, off_ijk, halo_width );
    }

    template <class DecompositionTag, std::size_t NSD = num_space_dim>
    std::enable_if_t<2 == NSD, IndexSpace<3>>
    sharedIndexSpace( DecompositionTag decomposition_tag, const int off_i,
                      const int off_j, const int halo_width = -1 ) const
    {
        std::array<int, 2> off_ijk = { off_i, off_j };
        return sharedIndexSpace( decomposition_tag, off_ijk, halo_width );
    }
 
  private:
    std::shared_ptr<LocalGrid<MeshType>> _local_grid;
    int _dofs_per_entity;
};

template <class>
struct is_array_layout : public std::false_type
{
};

template <class EntityType, class MeshType>
struct is_array_layout<ArrayLayout<EntityType, MeshType>>
    : public std::true_type
{
};

template <class EntityType, class MeshType>
struct is_array_layout<const ArrayLayout<EntityType, MeshType>>
    : public std::true_type
{
};
 
//---------------------------------------------------------------------------//
// Array layout creation.
//---------------------------------------------------------------------------//
template <class EntityType, class MeshType>
std::shared_ptr<ArrayLayout<EntityType, MeshType>>
createArrayLayout( const std::shared_ptr<LocalGrid<MeshType>>& local_grid,
                   const int dofs_per_entity, EntityType )
{
    return std::make_shared<ArrayLayout<EntityType, MeshType>>(
        local_grid, dofs_per_entity );
}
 
//---------------------------------------------------------------------------//
template <class EntityType, class MeshType>
std::shared_ptr<ArrayLayout<EntityType, MeshType>>
createArrayLayout( const std::shared_ptr<GlobalGrid<MeshType>>& global_grid,
                   const int halo_cell_width, const int dofs_per_entity,
                   EntityType )
{
    return std::make_shared<ArrayLayout<EntityType, MeshType>>(
        createLocalGrid( global_grid, halo_cell_width ), dofs_per_entity );
}
 
//---------------------------------------------------------------------------//
template <class Scalar, class EntityType, class MeshType, class... Params>
class Array
{
  public:
    using value_type = Scalar;

    using entity_type = EntityType;

    using mesh_type = MeshType;

    static constexpr std::size_t num_space_dim = mesh_type::num_space_dim;

    using array_layout = ArrayLayout<entity_type, mesh_type>;

    using view_type = std::conditional_t<
        3 == num_space_dim, Kokkos::View<value_type****, Params...>,
        std::conditional_t<2 == num_space_dim,
                           Kokkos::View<value_type***, Params...>, void>>;

    using memory_space = typename view_type::memory_space;
    using device_type [[deprecated]] = typename memory_space::device_type;
    using execution_space = typename memory_space::execution_space;

    Array( const std::string& label,
           const std::shared_ptr<array_layout>& layout )
        : _layout( layout )
        , _data( createView<value_type, Params...>(
              label, layout->indexSpace( Ghost(), Local() ) ) )
    {
    }

    Array( const std::shared_ptr<array_layout>& layout, const view_type& view )
        : _layout( layout )
        , _data( view )
    {
        for ( std::size_t d = 0; d < num_space_dim + 1; ++d )
            if ( (long)view.extent( d ) !=
                 layout->indexSpace( Ghost(), Local() ).extent( d ) )
                throw std::runtime_error(
                    "Cabana::Grid::Array: Layout and "
                    "view dimensions do not match (View label: " +
                    view.label() + ")" );
    }

    std::shared_ptr<array_layout> layout() const { return _layout; }

    view_type view() const { return _data; }

    std::string label() const { return _data.label(); }
 
  private:
    std::shared_ptr<array_layout> _layout;
    view_type _data;
 
  public:
    using subview_type = decltype( createSubview(
        _data, _layout->indexSpace( Ghost(), Local() ) ) );
    using subview_layout = typename subview_type::array_layout;
    using subview_memory_traits = typename subview_type::memory_traits;
    using subarray_type = Array<Scalar, EntityType, MeshType, subview_layout,
                                memory_space, subview_memory_traits>;
};
 
//---------------------------------------------------------------------------//
// Static type checker.
//---------------------------------------------------------------------------//
// Static type checker.
template <class>
struct is_array : public std::false_type
{
};
 
template <class Scalar, class EntityType, class MeshType, class... Params>
struct is_array<Array<Scalar, EntityType, MeshType, Params...>>
    : public std::true_type
{
};
 
template <class Scalar, class EntityType, class MeshType, class... Params>
struct is_array<const Array<Scalar, EntityType, MeshType, Params...>>
    : public std::true_type
{
};
 
//---------------------------------------------------------------------------//
// Array creation.
//---------------------------------------------------------------------------//
template <class Scalar, class... Params, class EntityType, class MeshType>
std::shared_ptr<Array<Scalar, EntityType, MeshType, Params...>>
createArray( const std::string& label,
             const std::shared_ptr<ArrayLayout<EntityType, MeshType>>& layout )
{
    return std::make_shared<Array<Scalar, EntityType, MeshType, Params...>>(
        label, layout );
}
 
//---------------------------------------------------------------------------//
// Subarray creation.
//---------------------------------------------------------------------------//
template <class Scalar, class EntityType, class MeshType, class... Params>
std::shared_ptr<Array<
    Scalar, EntityType, MeshType,
    typename Array<Scalar, EntityType, MeshType, Params...>::subview_layout,
    typename Array<Scalar, EntityType, MeshType, Params...>::memory_space,
    typename Array<Scalar, EntityType, MeshType,
                   Params...>::subview_memory_traits>>
createSubarray( const Array<Scalar, EntityType, MeshType, Params...>& array,
                const int dof_min, const int dof_max )
{
    if ( dof_min < 0 || dof_max > array.layout()->dofsPerEntity() )
        throw std::logic_error( "Cabana::Grid::createSubarray: Subarray "
                                "dimensions out of bounds (Label: " +
                                array.label() + ")" );
 
    auto space = array.layout()->indexSpace( Ghost(), Local() );
    std::array<long, MeshType::num_space_dim + 1> min;
    std::array<long, MeshType::num_space_dim + 1> max;
    for ( std::size_t d = 0; d < MeshType::num_space_dim; ++d )
    {
        min[d] = space.min( d );
        max[d] = space.max( d );
    }
    min.back() = dof_min;
    max.back() = dof_max;
    IndexSpace<MeshType::num_space_dim + 1> sub_space( min, max );
    auto sub_view = createSubview( array.view(), sub_space );
    auto sub_layout = createArrayLayout( array.layout()->localGrid(),
                                         dof_max - dof_min, EntityType() );
    return std::make_shared<Array<
        Scalar, EntityType, MeshType,
        typename Array<Scalar, EntityType, MeshType, Params...>::subview_layout,
        typename Array<Scalar, EntityType, MeshType, Params...>::memory_space,
        typename Array<Scalar, EntityType, MeshType,
                       Params...>::subview_memory_traits>>( sub_layout,
                                                            sub_view );
}
 
//---------------------------------------------------------------------------//
// Array operations.
//---------------------------------------------------------------------------//
namespace ArrayOp
{
//---------------------------------------------------------------------------//
template <class Scalar, class... Params, class EntityType, class MeshType>
std::shared_ptr<Array<Scalar, EntityType, MeshType, Params...>>
clone( const Array<Scalar, EntityType, MeshType, Params...>& array )
{
    return createArray<Scalar, Params...>( array.label(), array.layout() );
}
 
//---------------------------------------------------------------------------//
template <class Array_t, class DecompositionTag>
void assign( Array_t& array, const typename Array_t::value_type alpha,
             DecompositionTag tag )
{
    static_assert(
        is_array<Array_t>::value,
        "Cabana::Grid::ArrayOp::assign: Cabana::Grid::Array required" );
    auto subview = createSubview( array.view(),
                                  array.layout()->indexSpace( tag, Local() ) );
    Kokkos::deep_copy( subview, alpha );
}
 
//---------------------------------------------------------------------------//
template <class Array_t, class DecompositionTag>
std::enable_if_t<3 == Array_t::num_space_dim, void>
scale( Array_t& array, const typename Array_t::value_type alpha,
       DecompositionTag tag )
{
    static_assert(
        is_array<Array_t>::value,
        "Cabana::Grid::ArrayOp::scale: Cabana::Grid::Array required" );
    auto view = array.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::scale",
        createExecutionPolicy( array.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const int i, const int j, const int k, const int l ) {
            view( i, j, k, l ) *= alpha;
        } );
}

template <class Array_t, class DecompositionTag>
std::enable_if_t<2 == Array_t::num_space_dim, void>
scale( Array_t& array, const typename Array_t::value_type alpha,
       DecompositionTag tag )
{
    static_assert(
        is_array<Array_t>::value,
        "Cabana::Grid::ArrayOp::scale: Cabana::Grid::Array required" );
    auto view = array.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::scale",
        createExecutionPolicy( array.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const int i, const int j, const int l ) {
            view( i, j, l ) *= alpha;
        } );
}
 
//---------------------------------------------------------------------------//
template <class Array_t, class DecompositionTag>
std::enable_if_t<3 == Array_t::num_space_dim, void>
scale( Array_t& array, const std::vector<typename Array_t::value_type>& alpha,
       DecompositionTag tag )
{
    static_assert(
        is_array<Array_t>::value,
        "Cabana::Grid::ArrayOp::scale: Cabana::Grid::Array required" );
    if ( alpha.size() !=
         static_cast<unsigned>( array.layout()->dofsPerEntity() ) )
        throw std::runtime_error( "Cabana::Grid::ArrayOp::scale: Incorrect "
                                  "vector size (Label: " +
                                  array.label() + ")" );
 
    Kokkos::View<const typename Array_t::value_type*, Kokkos::HostSpace,
                 Kokkos::MemoryUnmanaged>
        alpha_view_host( alpha.data(), alpha.size() );
    auto alpha_view = Kokkos::create_mirror_view_and_copy(
        typename Array_t::memory_space(), alpha_view_host );
 
    auto array_view = array.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::scale",
        createExecutionPolicy( array.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const int i, const int j, const int k, const int l ) {
            array_view( i, j, k, l ) *= alpha_view( l );
        } );
}

template <class Array_t, class DecompositionTag>
std::enable_if_t<2 == Array_t::num_space_dim, void>
scale( Array_t& array, const std::vector<typename Array_t::value_type>& alpha,
       DecompositionTag tag )
{
    static_assert(
        is_array<Array_t>::value,
        "Cabana::Grid::ArrayOp::scale: Cabana::Grid::Array required" );
    if ( alpha.size() !=
         static_cast<unsigned>( array.layout()->dofsPerEntity() ) )
        throw std::runtime_error( "Cabana::Grid::ArrayOp::scale: Incorrect "
                                  "vector size (Label: " +
                                  array.label() + ")" );
 
    Kokkos::View<const typename Array_t::value_type*, Kokkos::HostSpace,
                 Kokkos::MemoryUnmanaged>
        alpha_view_host( alpha.data(), alpha.size() );
    auto alpha_view = Kokkos::create_mirror_view_and_copy(
        typename Array_t::memory_space(), alpha_view_host );
 
    auto array_view = array.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::scale",
        createExecutionPolicy( array.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const int i, const int j, const int l ) {
            array_view( i, j, l ) *= alpha_view( l );
        } );
}
 
//---------------------------------------------------------------------------//
template <class Array_t, class DecompositionTag>
void copy( Array_t& a, const Array_t& b, DecompositionTag tag )
{
    static_assert(
        is_array<Array_t>::value,
        "Cabana::Grid::ArrayOp::copy: Cabana::Grid::Array required" );
    auto a_space = a.layout()->indexSpace( tag, Local() );
    auto b_space = b.layout()->indexSpace( tag, Local() );
    if ( a_space != b_space )
        throw std::logic_error( "Cabana::Grid::ArrayOp::copy: Incompatible "
                                "index spaces (Labels: " +
                                a.label() + ", " + b.label() + ")" );
    auto subview_a = createSubview( a.view(), a_space );
    auto subview_b = createSubview( b.view(), b_space );
    Kokkos::deep_copy( subview_a, subview_b );
}
 
//---------------------------------------------------------------------------//
template <class Array_t, class DecompositionTag>
std::shared_ptr<Array_t> cloneCopy( const Array_t& array, DecompositionTag tag )
{
    auto cln = clone( array );
    copy( *cln, array, tag );
    return cln;
}
 
//---------------------------------------------------------------------------//
template <class Array_t, class DecompositionTag>
std::enable_if_t<3 == Array_t::num_space_dim, void>
update( Array_t& a, const typename Array_t::value_type alpha, const Array_t& b,
        const typename Array_t::value_type beta, DecompositionTag tag )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    auto a_view = a.view();
    auto b_view = b.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::update",
        createExecutionPolicy( a.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const int i, const int j, const int k, const int l ) {
            a_view( i, j, k, l ) =
                alpha * a_view( i, j, k, l ) + beta * b_view( i, j, k, l );
        } );
}

template <class Array_t, class DecompositionTag>
std::enable_if_t<2 == Array_t::num_space_dim, void>
update( Array_t& a, const typename Array_t::value_type alpha, const Array_t& b,
        const typename Array_t::value_type beta, DecompositionTag tag )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    auto a_view = a.view();
    auto b_view = b.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::update",
        createExecutionPolicy( a.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const long i, const long j, const long l ) {
            a_view( i, j, l ) =
                alpha * a_view( i, j, l ) + beta * b_view( i, j, l );
        } );
}
 
//---------------------------------------------------------------------------//
template <class Array_t, class DecompositionTag>
std::enable_if_t<3 == Array_t::num_space_dim, void>
update( Array_t& a, const typename Array_t::value_type alpha, const Array_t& b,
        const typename Array_t::value_type beta, const Array_t& c,
        const typename Array_t::value_type gamma, DecompositionTag tag )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    auto a_view = a.view();
    auto b_view = b.view();
    auto c_view = c.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::update",
        createExecutionPolicy( a.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const int i, const int j, const int k, const int l ) {
            a_view( i, j, k, l ) = alpha * a_view( i, j, k, l ) +
                                   beta * b_view( i, j, k, l ) +
                                   gamma * c_view( i, j, k, l );
        } );
}

template <class Array_t, class DecompositionTag>
std::enable_if_t<2 == Array_t::num_space_dim, void>
update( Array_t& a, const typename Array_t::value_type alpha, const Array_t& b,
        const typename Array_t::value_type beta, const Array_t& c,
        const typename Array_t::value_type gamma, DecompositionTag tag )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    auto a_view = a.view();
    auto b_view = b.view();
    auto c_view = c.view();
    Kokkos::parallel_for(
        "Cabana::Grid::ArrayOp::update",
        createExecutionPolicy( a.layout()->indexSpace( tag, Local() ),
                               typename Array_t::execution_space() ),
        KOKKOS_LAMBDA( const int i, const int j, const int l ) {
            a_view( i, j, l ) = alpha * a_view( i, j, l ) +
                                beta * b_view( i, j, l ) +
                                gamma * c_view( i, j, l );
        } );
}
 
//---------------------------------------------------------------------------//
template <class ViewType, std::size_t NumSpaceDim>
struct DotFunctor
{
    static constexpr std::size_t num_space_dim = NumSpaceDim;
    typedef typename ViewType::value_type value_type[];
    typedef typename ViewType::size_type size_type;
    size_type value_count;
    ViewType _a;
    ViewType _b;

    DotFunctor( const ViewType& a, const ViewType& b )
        : value_count( a.extent( NumSpaceDim ) )
        , _a( a )
        , _b( b )
    {
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<3 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type k,
                const size_type l, value_type sum ) const
    {
        sum[l] += _a( i, j, k, l ) * _b( i, j, k, l );
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<2 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type l,
                value_type sum ) const
    {
        sum[l] += _a( i, j, l ) * _b( i, j, l );
    }

    KOKKOS_INLINE_FUNCTION
    void join( value_type dst, const value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            dst[j] += src[j];
    }

    KOKKOS_INLINE_FUNCTION
    void join( volatile value_type dst, const volatile value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            dst[j] += src[j];
    }

    KOKKOS_INLINE_FUNCTION void init( value_type sum ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            sum[j] = 0.0;
    }
};

template <class Array_t>
void dot( const Array_t& a, const Array_t& b,
          std::vector<typename Array_t::value_type>& products )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    if ( products.size() !=
         static_cast<unsigned>( a.layout()->dofsPerEntity() ) )
        throw std::runtime_error(
            "Cabana::Grid::ArrayOp::dot: Incorrect vector size" );
 
    for ( auto& p : products )
        p = 0.0;
 
    DotFunctor<typename Array_t::view_type, Array_t::num_space_dim> functor(
        a.view(), b.view() );
    typename Array_t::execution_space exec_space;
    Kokkos::parallel_reduce(
        "Cabana::Grid::ArrayOp::dot",
        createExecutionPolicy( a.layout()->indexSpace( Own(), Local() ),
                               exec_space ),
        functor,
        Kokkos::View<typename Array_t::value_type*, Kokkos::HostSpace>(
            products.data(), products.size() ) );
    exec_space.fence( "Cabana::Grid::ArrayOp::dot before MPI_Allreduce" );
 
    MPI_Allreduce( MPI_IN_PLACE, products.data(), products.size(),
                   MpiTraits<typename Array_t::value_type>::type(), MPI_SUM,
                   a.layout()->localGrid()->globalGrid().comm() );
}
 
//---------------------------------------------------------------------------//
template <class ViewType, std::size_t NumSpaceDim>
struct NormInfFunctor
{
    static constexpr std::size_t num_space_dim = NumSpaceDim;
    typedef typename ViewType::value_type value_type[];
    typedef typename ViewType::size_type size_type;
    size_type value_count;
    ViewType _view;

    NormInfFunctor( const ViewType& view )
        : value_count( view.extent( NumSpaceDim ) )
        , _view( view )
    {
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<3 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type k,
                const size_type l, value_type norm ) const
    {
        auto v_abs = fabs( _view( i, j, k, l ) );
        if ( v_abs > norm[l] )
            norm[l] = v_abs;
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<2 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type l,
                value_type norm ) const
    {
        auto v_abs = fabs( _view( i, j, l ) );
        if ( v_abs > norm[l] )
            norm[l] = v_abs;
    }

    KOKKOS_INLINE_FUNCTION
    void join( value_type dst, const value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            if ( src[j] > dst[j] )
                dst[j] = src[j];
    }

    KOKKOS_INLINE_FUNCTION
    void join( volatile value_type dst, const volatile value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            if ( src[j] > dst[j] )
                dst[j] = src[j];
    }

    KOKKOS_INLINE_FUNCTION void init( value_type norm ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            norm[j] = 0.0;
    }
};

template <class Array_t>
void normInf( const Array_t& array,
              std::vector<typename Array_t::value_type>& norms )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    if ( norms.size() !=
         static_cast<unsigned>( array.layout()->dofsPerEntity() ) )
        throw std::runtime_error(
            "Cabana::Grid::ArrayOp::normInf: Incorrect vector size" );
 
    for ( auto& n : norms )
        n = 0.0;
 
    NormInfFunctor<typename Array_t::view_type, Array_t::num_space_dim> functor(
        array.view() );
    typename Array_t::execution_space exec_space;
    Kokkos::parallel_reduce(
        "Cabana::Grid::ArrayOp::normInf",
        createExecutionPolicy( array.layout()->indexSpace( Own(), Local() ),
                               exec_space ),
        functor,
        Kokkos::View<typename Array_t::value_type*, Kokkos::HostSpace>(
            norms.data(), norms.size() ) );
    exec_space.fence( "ArrayOp::normInf before MPI_Allreduce" );
 
    MPI_Allreduce( MPI_IN_PLACE, norms.data(), norms.size(),
                   MpiTraits<typename Array_t::value_type>::type(), MPI_MAX,
                   array.layout()->localGrid()->globalGrid().comm() );
}
 
//---------------------------------------------------------------------------//
template <class ViewType, std::size_t NumSpaceDim>
struct Norm1Functor
{
    static constexpr std::size_t num_space_dim = NumSpaceDim;
    typedef typename ViewType::value_type value_type[];
    typedef typename ViewType::size_type size_type;
    size_type value_count;
    ViewType _view;

    Norm1Functor( const ViewType& view )
        : value_count( view.extent( NumSpaceDim ) )
        , _view( view )
    {
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<3 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type k,
                const size_type l, value_type norm ) const
    {
        norm[l] += fabs( _view( i, j, k, l ) );
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<2 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type l,
                value_type norm ) const
    {
        norm[l] += fabs( _view( i, j, l ) );
    }

    KOKKOS_INLINE_FUNCTION
    void join( value_type dst, const value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            dst[j] += src[j];
    }

    KOKKOS_INLINE_FUNCTION
    void join( volatile value_type dst, const volatile value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            dst[j] += src[j];
    }

    KOKKOS_INLINE_FUNCTION void init( value_type norm ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            norm[j] = 0.0;
    }
};

template <class Array_t>
void norm1( const Array_t& array,
            std::vector<typename Array_t::value_type>& norms )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    if ( norms.size() !=
         static_cast<unsigned>( array.layout()->dofsPerEntity() ) )
        throw std::runtime_error(
            "Cabana::Grid::ArrayOp::norm1: Incorrect vector size" );
 
    for ( auto& n : norms )
        n = 0.0;
 
    Norm1Functor<typename Array_t::view_type, Array_t::num_space_dim> functor(
        array.view() );
    typename Array_t::execution_space exec_space;
    Kokkos::parallel_reduce(
        "Cabana::Grid::ArrayOp::norm1",
        createExecutionPolicy( array.layout()->indexSpace( Own(), Local() ),
                               exec_space ),
        functor,
        Kokkos::View<typename Array_t::value_type*, Kokkos::HostSpace>(
            norms.data(), norms.size() ) );
    exec_space.fence( "Cabana::Grid::ArrayOp::norm1 before MPI_Allreduce" );
 
    MPI_Allreduce( MPI_IN_PLACE, norms.data(), norms.size(),
                   MpiTraits<typename Array_t::value_type>::type(), MPI_SUM,
                   array.layout()->localGrid()->globalGrid().comm() );
}
 
//---------------------------------------------------------------------------//
template <class ViewType, std::size_t NumSpaceDim>
struct Norm2Functor
{
    static constexpr std::size_t num_space_dim = NumSpaceDim;
    typedef typename ViewType::value_type value_type[];
    typedef typename ViewType::size_type size_type;
    size_type value_count;
    ViewType _view;

    Norm2Functor( const ViewType& view )
        : value_count( view.extent( NumSpaceDim ) )
        , _view( view )
    {
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<3 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type k,
                const size_type l, value_type norm ) const
    {
        norm[l] += _view( i, j, k, l ) * _view( i, j, k, l );
    }

    template <std::size_t NSD = num_space_dim>
    KOKKOS_INLINE_FUNCTION std::enable_if_t<2 == NSD, void>
    operator()( const size_type i, const size_type j, const size_type l,
                value_type norm ) const
    {
        norm[l] += _view( i, j, l ) * _view( i, j, l );
    }

    KOKKOS_INLINE_FUNCTION
    void join( value_type dst, const value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            dst[j] += src[j];
    }

    KOKKOS_INLINE_FUNCTION
    void join( volatile value_type dst, const volatile value_type src ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            dst[j] += src[j];
    }

    KOKKOS_INLINE_FUNCTION void init( value_type norm ) const
    {
        for ( size_type j = 0; j < value_count; ++j )
            norm[j] = 0.0;
    }
};

template <class Array_t>
void norm2( const Array_t& array,
            std::vector<typename Array_t::value_type>& norms )
{
    static_assert( is_array<Array_t>::value, "Cabana::Grid::Array required" );
    if ( norms.size() !=
         static_cast<unsigned>( array.layout()->dofsPerEntity() ) )
        throw std::runtime_error(
            "Cabana::Grid::ArrayOp::norm2: Incorrect vector size" );
 
    for ( auto& n : norms )
        n = 0.0;
 
    Norm2Functor<typename Array_t::view_type, Array_t::num_space_dim> functor(
        array.view() );
    typename Array_t::execution_space exec_space;
    Kokkos::parallel_reduce(
        "Cabana::Grid::ArrayOp::norm2",
        createExecutionPolicy( array.layout()->indexSpace( Own(), Local() ),
                               exec_space ),
        functor,
        Kokkos::View<typename Array_t::value_type*, Kokkos::HostSpace>(
            norms.data(), norms.size() ) );
    exec_space.fence( "Cabana::Grid::ArrayOp::norm2 before MPI_Allreduce" );
 
    MPI_Allreduce( MPI_IN_PLACE, norms.data(), norms.size(),
                   MpiTraits<typename Array_t::value_type>::type(), MPI_SUM,
                   array.layout()->localGrid()->globalGrid().comm() );
 
    for ( auto& n : norms )
        n = std::sqrt( n );
}
//---------------------------------------------------------------------------//
 
} // end namespace ArrayOp
 
//---------------------------------------------------------------------------//
 
} // namespace Grid
} // namespace Cabana
 
#endif // end CABANA_GRID_ARRAY_HPP