From d9cc46855c2f254aa1a2c7fe67c6787e0ffad0a9 Mon Sep 17 00:00:00 2001 From: Sam Reeve <6740307+streeve@users.noreply.github.com> Date: Fri, 17 Feb 2023 17:29:32 -0500 Subject: [PATCH] Test both Views and slices for interpolation --- cajita/unit_test/tstInterpolation2d.hpp | 444 +++++++++++++++++------ cajita/unit_test/tstInterpolation3d.hpp | 450 ++++++++++++++++++------ 2 files changed, 665 insertions(+), 229 deletions(-) diff --git a/cajita/unit_test/tstInterpolation2d.hpp b/cajita/unit_test/tstInterpolation2d.hpp index 480705c72..e82c1d61e 100644 --- a/cajita/unit_test/tstInterpolation2d.hpp +++ b/cajita/unit_test/tstInterpolation2d.hpp @@ -11,6 +11,10 @@ #include +#include +#include +#include + #include #include #include @@ -34,97 +38,22 @@ using namespace Cajita; namespace Test { -//---------------------------------------------------------------------------// -void interpolationTest() +template +void checkP2G( const PositionType points, ScalarPoint scalar_point_field, + VectorPoint vector_point_field, TensorPoint tensor_point_field, + ScalarGrid scalar_grid_field, VectorGrid vector_grid_field, + const NodeIndexSpace node_space, const ScalarHalo scalar_halo, + const VectorHalo vector_halo, const std::size_t num_point ) { - // Create the global mesh. - std::array low_corner = { -1.2, 0.1 }; - std::array high_corner = { -0.2, 9.5 }; - double cell_size = 0.1; - auto global_mesh = - createUniformGlobalMesh( low_corner, high_corner, cell_size ); - - // Create the global grid. - DimBlockPartitioner<2> partitioner; - std::array is_dim_periodic = { true, true }; - auto global_grid = createGlobalGrid( MPI_COMM_WORLD, global_mesh, - is_dim_periodic, partitioner ); - - // Create a grid local_grid. - int halo_width = 1; - auto local_grid = createLocalGrid( global_grid, halo_width ); - auto local_mesh = createLocalMesh( *local_grid ); - - // Create a point in the center of every cell. - auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); - int num_point = cell_space.size(); - Kokkos::View points( - Kokkos::ViewAllocateWithoutInitializing( "points" ), num_point ); - Kokkos::parallel_for( - "fill_points", createExecutionPolicy( cell_space, TEST_EXECSPACE() ), - KOKKOS_LAMBDA( const int i, const int j ) { - int pi = i - halo_width; - int pj = j - halo_width; - int pid = pi + cell_space.extent( Dim::I ) * pj; - int idx[2] = { i, j }; - double x[2]; - local_mesh.coordinates( Cell(), idx, x ); - points( pid, Dim::I ) = x[Dim::I]; - points( pid, Dim::J ) = x[Dim::J]; - } ); - - // Node space. - auto node_space = local_grid->indexSpace( Own(), Node(), Local() ); - - // Create a scalar field on the grid. - auto scalar_layout = createArrayLayout( local_grid, 1, Node() ); - auto scalar_grid_field = - createArray( "scalar_grid_field", scalar_layout ); - auto scalar_halo = - createHalo( NodeHaloPattern<2>(), halo_width, *scalar_grid_field ); - auto scalar_grid_host = - Kokkos::create_mirror_view( scalar_grid_field->view() ); - - // Create a vector field on the grid. - auto vector_layout = createArrayLayout( local_grid, 2, Node() ); - auto vector_grid_field = - createArray( "vector_grid_field", vector_layout ); - auto vector_halo = - createHalo( NodeHaloPattern<2>(), halo_width, *vector_grid_field ); - auto vector_grid_host = - Kokkos::create_mirror_view( vector_grid_field->view() ); - - // Create a scalar point field. - Kokkos::View scalar_point_field( - Kokkos::ViewAllocateWithoutInitializing( "scalar_point_field" ), - num_point ); - auto scalar_point_host = Kokkos::create_mirror_view( scalar_point_field ); - - // Create a vector point field. - Kokkos::View vector_point_field( - Kokkos::ViewAllocateWithoutInitializing( "vector_point_field" ), - num_point ); - auto vector_point_host = Kokkos::create_mirror_view( vector_point_field ); - - // Create a tensor point field. - Kokkos::View tensor_point_field( - Kokkos::ViewAllocateWithoutInitializing( "tensor_point_field" ), - num_point ); - auto tensor_point_host = Kokkos::create_mirror_view( tensor_point_field ); - - // P2G - // --- - - Kokkos::deep_copy( scalar_point_field, 3.5 ); - Kokkos::deep_copy( vector_point_field, 3.5 ); - Kokkos::deep_copy( tensor_point_field, 3.5 ); - // Interpolate a scalar point value to the grid. ArrayOp::assign( *scalar_grid_field, 0.0, Ghost() ); auto scalar_p2g = createScalarValueP2G( scalar_point_field, -0.5 ); p2g( scalar_p2g, points, num_point, Spline<1>(), *scalar_halo, *scalar_grid_field ); - Kokkos::deep_copy( scalar_grid_host, scalar_grid_field->view() ); + auto scalar_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -135,7 +64,8 @@ void interpolationTest() auto vector_p2g = createVectorValueP2G( vector_point_field, -0.5 ); p2g( vector_p2g, points, num_point, Spline<1>(), *vector_halo, *vector_grid_field ); - Kokkos::deep_copy( vector_grid_host, vector_grid_field->view() ); + auto vector_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -147,7 +77,8 @@ void interpolationTest() auto scalar_grad_p2g = createScalarGradientP2G( scalar_point_field, -0.5 ); p2g( scalar_grad_p2g, points, num_point, Spline<1>(), *vector_halo, *vector_grid_field ); - Kokkos::deep_copy( vector_grid_host, vector_grid_field->view() ); + vector_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -159,7 +90,8 @@ void interpolationTest() auto vector_div_p2g = createVectorDivergenceP2G( vector_point_field, -0.5 ); p2g( vector_div_p2g, points, num_point, Spline<1>(), *scalar_halo, *scalar_grid_field ); - Kokkos::deep_copy( scalar_grid_host, scalar_grid_field->view() ); + scalar_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -170,75 +102,361 @@ void interpolationTest() auto tensor_div_p2g = createTensorDivergenceP2G( tensor_point_field, -0.5 ); p2g( tensor_div_p2g, points, num_point, Spline<1>(), *vector_halo, *vector_grid_field ); - Kokkos::deep_copy( vector_grid_host, vector_grid_field->view() ); + vector_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) for ( int d = 0; d < 2; ++d ) EXPECT_FLOAT_EQ( vector_grid_host( i, j, d ) + 1.0, 1.0 ); +} - // G2P - // --- - +template +void scalarValueG2P( const PositionType points, ScalarPoint& scalar_point_field, + ScalarGrid& scalar_grid_field, + const ScalarHalo scalar_halo, const std::size_t num_point ) +{ ArrayOp::assign( *scalar_grid_field, 3.5, Own() ); - ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a scalar grid value to the points. - Kokkos::deep_copy( scalar_point_field, 0.0 ); auto scalar_value_g2p = createScalarValueG2P( scalar_point_field, -0.5 ); g2p( *scalar_grid_field, *scalar_halo, points, num_point, Spline<1>(), scalar_value_g2p ); - Kokkos::deep_copy( scalar_point_host, scalar_point_field ); - for ( int p = 0; p < num_point; ++p ) +} + +template +void checkScalarValueG2P( const ScalarPoint scalar_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) EXPECT_FLOAT_EQ( scalar_point_host( p ), -1.75 ); +} + +template +void vectorValueG2P( const PositionType points, VectorPoint& vector_point_field, + VectorGrid vector_grid_field, const VectorHalo vector_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a vector grid value to the points. - Kokkos::deep_copy( vector_point_field, 0.0 ); auto vector_value_g2p = createVectorValueG2P( vector_point_field, -0.5 ); g2p( *vector_grid_field, *vector_halo, points, num_point, Spline<1>(), vector_value_g2p ); - Kokkos::deep_copy( vector_point_host, vector_point_field ); - for ( int p = 0; p < num_point; ++p ) - for ( int d = 0; d < 2; ++d ) +} + +template +void checkVectorValueG2P( const VectorPoint vector_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) + for ( std::size_t d = 0; d < 2; ++d ) EXPECT_FLOAT_EQ( vector_point_host( p, d ), -1.75 ); +} + +template +void scalarGradientG2P( const PositionType points, + VectorPoint& vector_point_field, + ScalarGrid scalar_grid_field, + const ScalarHalo scalar_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *scalar_grid_field, 3.5, Own() ); // Interpolate a scalar grid gradient to the points. - Kokkos::deep_copy( vector_point_field, 0.0 ); auto scalar_gradient_g2p = createScalarGradientG2P( vector_point_field, -0.5 ); g2p( *scalar_grid_field, *scalar_halo, points, num_point, Spline<1>(), scalar_gradient_g2p ); - Kokkos::deep_copy( vector_point_host, vector_point_field ); - for ( int p = 0; p < num_point; ++p ) - for ( int d = 0; d < 2; ++d ) +} + +template +void checkScalarGradientG2P( const VectorPoint vector_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) + for ( std::size_t d = 0; d < 2; ++d ) EXPECT_FLOAT_EQ( vector_point_host( p, d ) + 1.0, 1.0 ); +} + +template +void vectorGradientG2P( const PositionType points, + TensorPoint& tensor_point_field, + VectorGrid vector_grid_field, + const VectorHalo vector_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a vector grid gradient to the points. - Kokkos::deep_copy( tensor_point_field, 0.0 ); auto vector_gradient_g2p = createVectorGradientG2P( tensor_point_field, -0.5 ); g2p( *vector_grid_field, *vector_halo, points, num_point, Spline<1>(), vector_gradient_g2p ); - Kokkos::deep_copy( tensor_point_host, tensor_point_field ); - for ( int p = 0; p < num_point; ++p ) - for ( int i = 0; i < 2; ++i ) - for ( int j = 0; j < 2; ++j ) +} + +template +void checkVectorGradientG2P( const TensorPoint tensor_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) + for ( std::size_t i = 0; i < 2; ++i ) + for ( std::size_t j = 0; j < 2; ++j ) EXPECT_FLOAT_EQ( tensor_point_host( p, i, j ) + 1.0, 1.0 ); +} + +template +void vectorDivergenceG2P( const PositionType points, + ScalarPoint& scalar_point_field, + VectorGrid vector_grid_field, + const VectorHalo vector_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a vector grid divergence to the points. - Kokkos::deep_copy( scalar_point_field, 0.0 ); auto vector_div_g2p = createVectorDivergenceG2P( scalar_point_field, -0.5 ); g2p( *vector_grid_field, *vector_halo, points, num_point, Spline<1>(), vector_div_g2p ); - Kokkos::deep_copy( scalar_point_host, scalar_point_field ); - for ( int p = 0; p < num_point; ++p ) +} + +template +void checkVectorDivergenceG2P( const ScalarPoint scalar_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) EXPECT_FLOAT_EQ( scalar_point_host( p ) + 1.0, 1.0 ); } +auto createGrid() +{ + // Create the global mesh. + std::array low_corner = { -1.2, 0.1 }; + std::array high_corner = { -0.2, 9.5 }; + double cell_size = 0.1; + auto global_mesh = + createUniformGlobalMesh( low_corner, high_corner, cell_size ); + + // Create the global grid. + DimBlockPartitioner<2> partitioner; + std::array is_dim_periodic = { true, true }; + auto global_grid = createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_dim_periodic, partitioner ); + return global_grid; +} + +template +void setPositions( PositionType& points, const LocalGridType local_grid, + const int halo_width ) +{ + auto local_mesh = createLocalMesh( *local_grid ); + auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); + + // Create a point in the center of every cell. + Kokkos::parallel_for( + "fill_points", createExecutionPolicy( cell_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j ) { + int pi = i - halo_width; + int pj = j - halo_width; + int pid = pi + cell_space.extent( Dim::I ) * pj; + int idx[2] = { i, j }; + double x[2]; + local_mesh.coordinates( Cell(), idx, x ); + points( pid, Dim::I ) = x[Dim::I]; + points( pid, Dim::J ) = x[Dim::J]; + } ); +} +//---------------------------------------------------------------------------// +void interpolationViewTest() +{ + // Create a grid local_grid. + int halo_width = 1; + auto global_grid = createGrid(); + auto local_grid = createLocalGrid( global_grid, halo_width ); + auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); + auto num_point = cell_space.size(); + + // Create particles. + Kokkos::View points( + Kokkos::ViewAllocateWithoutInitializing( "points" ), num_point ); + setPositions( points, local_grid, halo_width ); + + // Node space. + auto node_space = local_grid->indexSpace( Own(), Node(), Local() ); + + // Create a scalar field on the grid. + auto scalar_layout = createArrayLayout( local_grid, 1, Node() ); + auto scalar_grid_field = + createArray( "scalar_grid_field", scalar_layout ); + + // Create a vector field on the grid. + auto vector_layout = createArrayLayout( local_grid, 2, Node() ); + auto vector_grid_field = + createArray( "vector_grid_field", vector_layout ); + + // Create a scalar, vector, and tensor point fields. + Kokkos::View scalar_point_field( + Kokkos::ViewAllocateWithoutInitializing( "scalar_point_field" ), + num_point ); + Kokkos::View vector_point_field( + Kokkos::ViewAllocateWithoutInitializing( "vector_point_field" ), + num_point ); + Kokkos::View tensor_point_field( + Kokkos::ViewAllocateWithoutInitializing( "tensor_point_field" ), + num_point ); + + // Create halos. + auto scalar_halo = + createHalo( NodeHaloPattern<2>(), halo_width, *scalar_grid_field ); + auto vector_halo = + createHalo( NodeHaloPattern<2>(), halo_width, *vector_grid_field ); + + Kokkos::deep_copy( scalar_point_field, 3.5 ); + Kokkos::deep_copy( vector_point_field, 3.5 ); + Kokkos::deep_copy( tensor_point_field, 3.5 ); + + // Check all p2g together. + checkP2G( points, scalar_point_field, vector_point_field, + tensor_point_field, scalar_grid_field, vector_grid_field, + node_space, scalar_halo, vector_halo, num_point ); + + // Separated G2P because of reset/copies (different for AoSoA/View). + Kokkos::deep_copy( scalar_point_field, 0.0 ); + scalarValueG2P( points, scalar_point_field, scalar_grid_field, scalar_halo, + num_point ); + auto scalar_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_point_field ); + checkScalarValueG2P( scalar_point_host, num_point ); + + Kokkos::deep_copy( vector_point_field, 0.0 ); + vectorValueG2P( points, vector_point_field, vector_grid_field, vector_halo, + num_point ); + auto vector_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_point_field ); + checkVectorValueG2P( vector_point_host, num_point ); + + Kokkos::deep_copy( vector_point_field, 0.0 ); + scalarGradientG2P( points, vector_point_field, scalar_grid_field, + scalar_halo, num_point ); + vector_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_point_field ); + checkScalarGradientG2P( vector_point_host, num_point ); + + Kokkos::deep_copy( tensor_point_field, 0.0 ); + vectorGradientG2P( points, tensor_point_field, vector_grid_field, + vector_halo, num_point ); + auto tensor_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), tensor_point_field ); + checkVectorGradientG2P( tensor_point_host, num_point ); + + Kokkos::deep_copy( scalar_point_field, 0.0 ); + vectorDivergenceG2P( points, scalar_point_field, vector_grid_field, + vector_halo, num_point ); + scalar_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_point_field ); + checkVectorDivergenceG2P( scalar_point_host, num_point ); +} + +//---------------------------------------------------------------------------// +void interpolationSliceTest() +{ + // Create a grid local_grid. + int halo_width = 1; + auto global_grid = createGrid(); + auto local_grid = createLocalGrid( global_grid, halo_width ); + auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); + auto num_point = cell_space.size(); + + // Create particles. + using member_types = + Cabana::MemberTypes; + Cabana::AoSoA particles( "particles", + num_point ); + auto points = Cabana::slice<2>( particles, "position" ); + setPositions( points, local_grid, halo_width ); + + // Node space. + auto node_space = local_grid->indexSpace( Own(), Node(), Local() ); + + // Create a scalar field on the grid. + auto scalar_layout = createArrayLayout( local_grid, 1, Node() ); + auto scalar_grid_field = + createArray( "scalar_grid_field", scalar_layout ); + + // Create a vector field on the grid. + auto vector_layout = createArrayLayout( local_grid, 2, Node() ); + auto vector_grid_field = + createArray( "vector_grid_field", vector_layout ); + + // Create a scalar, vector, and tensor point fields. + auto tensor_point_field = Cabana::slice<0>( particles, "tensor" ); + auto vector_point_field = Cabana::slice<1>( particles, "vector" ); + auto scalar_point_field = Cabana::slice<3>( particles, "scalar" ); + + // Create halos. + auto scalar_halo = + createHalo( NodeHaloPattern<2>(), halo_width, *scalar_grid_field ); + auto vector_halo = + createHalo( NodeHaloPattern<2>(), halo_width, *vector_grid_field ); + + Cabana::deep_copy( scalar_point_field, 3.5 ); + Cabana::deep_copy( vector_point_field, 3.5 ); + Cabana::deep_copy( tensor_point_field, 3.5 ); + + // Check all p2g together. + checkP2G( points, scalar_point_field, vector_point_field, + tensor_point_field, scalar_grid_field, vector_grid_field, + node_space, scalar_halo, vector_halo, num_point ); + + // Separated to do the reset/copies (different for AoSoA/View). + auto particles_host = + Cabana::create_mirror_view( Kokkos::HostSpace(), particles ); + auto tensor_point_host = Cabana::slice<0>( particles, "tensor" ); + auto vector_point_host = Cabana::slice<1>( particles, "vector" ); + auto scalar_point_host = Cabana::slice<3>( particles, "scalar" ); + + Cabana::deep_copy( scalar_point_field, 0.0 ); + scalarValueG2P( points, scalar_point_field, scalar_grid_field, scalar_halo, + num_point ); + Cabana::deep_copy( scalar_point_host, scalar_point_field ); + checkScalarValueG2P( scalar_point_host, num_point ); + + Cabana::deep_copy( vector_point_field, 0.0 ); + vectorValueG2P( points, vector_point_field, vector_grid_field, vector_halo, + num_point ); + Cabana::deep_copy( vector_point_host, vector_point_field ); + checkVectorValueG2P( vector_point_host, num_point ); + + Cabana::deep_copy( vector_point_field, 0.0 ); + scalarGradientG2P( points, vector_point_field, scalar_grid_field, + scalar_halo, num_point ); + Cabana::deep_copy( vector_point_host, vector_point_field ); + checkScalarGradientG2P( vector_point_host, num_point ); + + Cabana::deep_copy( tensor_point_field, 0.0 ); + vectorGradientG2P( points, tensor_point_field, vector_grid_field, + vector_halo, num_point ); + Cabana::deep_copy( tensor_point_host, tensor_point_field ); + checkVectorGradientG2P( tensor_point_host, num_point ); + + Cabana::deep_copy( scalar_point_field, 0.0 ); + vectorDivergenceG2P( points, scalar_point_field, vector_grid_field, + vector_halo, num_point ); + Cabana::deep_copy( scalar_point_host, scalar_point_field ); + checkVectorDivergenceG2P( scalar_point_host, num_point ); +} + //---------------------------------------------------------------------------// // RUN TESTS //---------------------------------------------------------------------------// -TEST( interpolation, interpolation_test ) { interpolationTest(); } +TEST( interpolation, interpolation_view_test ) { interpolationViewTest(); } + +TEST( interpolation, interpolation_slice_test ) { interpolationSliceTest(); } //---------------------------------------------------------------------------// diff --git a/cajita/unit_test/tstInterpolation3d.hpp b/cajita/unit_test/tstInterpolation3d.hpp index 1e3e9a16f..83c18e19c 100644 --- a/cajita/unit_test/tstInterpolation3d.hpp +++ b/cajita/unit_test/tstInterpolation3d.hpp @@ -11,6 +11,10 @@ #include +#include +#include +#include + #include #include #include @@ -34,100 +38,22 @@ using namespace Cajita; namespace Test { -//---------------------------------------------------------------------------// -void interpolationTest() +template +void checkP2G( const PositionType points, ScalarPoint scalar_point_field, + VectorPoint vector_point_field, TensorPoint tensor_point_field, + ScalarGrid scalar_grid_field, VectorGrid vector_grid_field, + const NodeIndexSpace node_space, const ScalarHalo scalar_halo, + const VectorHalo vector_halo, const std::size_t num_point ) { - // Create the global mesh. - std::array low_corner = { -1.2, 0.1, 1.1 }; - std::array high_corner = { -0.3, 9.5, 2.3 }; - double cell_size = 0.1; - auto global_mesh = - createUniformGlobalMesh( low_corner, high_corner, cell_size ); - - // Create the global grid. - DimBlockPartitioner<3> partitioner; - std::array is_dim_periodic = { true, true, true }; - auto global_grid = createGlobalGrid( MPI_COMM_WORLD, global_mesh, - is_dim_periodic, partitioner ); - - // Create a grid local_grid. - int halo_width = 1; - auto local_grid = createLocalGrid( global_grid, halo_width ); - auto local_mesh = createLocalMesh( *local_grid ); - - // Create a point in the center of every cell. - auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); - int num_point = cell_space.size(); - Kokkos::View points( - Kokkos::ViewAllocateWithoutInitializing( "points" ), num_point ); - Kokkos::parallel_for( - "fill_points", createExecutionPolicy( cell_space, TEST_EXECSPACE() ), - KOKKOS_LAMBDA( const int i, const int j, const int k ) { - int pi = i - halo_width; - int pj = j - halo_width; - int pk = k - halo_width; - int pid = pi + cell_space.extent( Dim::I ) * - ( pj + cell_space.extent( Dim::J ) * pk ); - int idx[3] = { i, j, k }; - double x[3]; - local_mesh.coordinates( Cell(), idx, x ); - points( pid, Dim::I ) = x[Dim::I]; - points( pid, Dim::J ) = x[Dim::J]; - points( pid, Dim::K ) = x[Dim::K]; - } ); - - // Node space. - auto node_space = local_grid->indexSpace( Own(), Node(), Local() ); - - // Create a scalar field on the grid. - auto scalar_layout = createArrayLayout( local_grid, 1, Node() ); - auto scalar_grid_field = - createArray( "scalar_grid_field", scalar_layout ); - auto scalar_halo = - createHalo( NodeHaloPattern<3>(), halo_width, *scalar_grid_field ); - auto scalar_grid_host = - Kokkos::create_mirror_view( scalar_grid_field->view() ); - - // Create a vector field on the grid. - auto vector_layout = createArrayLayout( local_grid, 3, Node() ); - auto vector_grid_field = - createArray( "vector_grid_field", vector_layout ); - auto vector_halo = - createHalo( NodeHaloPattern<3>(), halo_width, *vector_grid_field ); - auto vector_grid_host = - Kokkos::create_mirror_view( vector_grid_field->view() ); - - // Create a scalar point field. - Kokkos::View scalar_point_field( - Kokkos::ViewAllocateWithoutInitializing( "scalar_point_field" ), - num_point ); - auto scalar_point_host = Kokkos::create_mirror_view( scalar_point_field ); - - // Create a vector point field. - Kokkos::View vector_point_field( - Kokkos::ViewAllocateWithoutInitializing( "vector_point_field" ), - num_point ); - auto vector_point_host = Kokkos::create_mirror_view( vector_point_field ); - - // Create a tensor point field. - Kokkos::View tensor_point_field( - Kokkos::ViewAllocateWithoutInitializing( "tensor_point_field" ), - num_point ); - auto tensor_point_host = Kokkos::create_mirror_view( tensor_point_field ); - - // P2G - // --- - - Kokkos::deep_copy( scalar_point_field, 3.5 ); - Kokkos::deep_copy( vector_point_field, 3.5 ); - Kokkos::deep_copy( tensor_point_field, 3.5 ); - // Interpolate a scalar point value to the grid. ArrayOp::assign( *scalar_grid_field, 0.0, Ghost() ); auto scalar_p2g = createScalarValueP2G( scalar_point_field, -0.5 ); p2g( scalar_p2g, points, num_point, Spline<1>(), *scalar_halo, *scalar_grid_field ); - Kokkos::deep_copy( scalar_grid_host, scalar_grid_field->view() ); + auto scalar_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -140,7 +66,8 @@ void interpolationTest() auto vector_p2g = createVectorValueP2G( vector_point_field, -0.5 ); p2g( vector_p2g, points, num_point, Spline<1>(), *vector_halo, *vector_grid_field ); - Kokkos::deep_copy( vector_grid_host, vector_grid_field->view() ); + auto vector_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -154,7 +81,8 @@ void interpolationTest() auto scalar_grad_p2g = createScalarGradientP2G( scalar_point_field, -0.5 ); p2g( scalar_grad_p2g, points, num_point, Spline<1>(), *vector_halo, *vector_grid_field ); - Kokkos::deep_copy( vector_grid_host, vector_grid_field->view() ); + vector_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -169,7 +97,8 @@ void interpolationTest() auto vector_div_p2g = createVectorDivergenceP2G( vector_point_field, -0.5 ); p2g( vector_div_p2g, points, num_point, Spline<1>(), *scalar_halo, *scalar_grid_field ); - Kokkos::deep_copy( scalar_grid_host, scalar_grid_field->view() ); + scalar_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -182,7 +111,8 @@ void interpolationTest() auto tensor_div_p2g = createTensorDivergenceP2G( tensor_point_field, -0.5 ); p2g( tensor_div_p2g, points, num_point, Spline<1>(), *vector_halo, *vector_grid_field ); - Kokkos::deep_copy( vector_grid_host, vector_grid_field->view() ); + vector_grid_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_grid_field->view() ); for ( int i = node_space.min( Dim::I ); i < node_space.max( Dim::I ); ++i ) for ( int j = node_space.min( Dim::J ); j < node_space.max( Dim::J ); ++j ) @@ -191,69 +121,357 @@ void interpolationTest() for ( int d = 0; d < 3; ++d ) EXPECT_FLOAT_EQ( vector_grid_host( i, j, k, d ) + 1.0, 1.0 ); +} - // G2P - // --- - +template +void scalarValueG2P( const PositionType points, ScalarPoint& scalar_point_field, + ScalarGrid& scalar_grid_field, + const ScalarHalo scalar_halo, const std::size_t num_point ) +{ ArrayOp::assign( *scalar_grid_field, 3.5, Own() ); - ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a scalar grid value to the points. - Kokkos::deep_copy( scalar_point_field, 0.0 ); auto scalar_value_g2p = createScalarValueG2P( scalar_point_field, -0.5 ); g2p( *scalar_grid_field, *scalar_halo, points, num_point, Spline<1>(), scalar_value_g2p ); - Kokkos::deep_copy( scalar_point_host, scalar_point_field ); - for ( int p = 0; p < num_point; ++p ) +} + +template +void checkScalarValueG2P( const ScalarPoint scalar_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) EXPECT_FLOAT_EQ( scalar_point_host( p ), -1.75 ); +} + +template +void vectorValueG2P( const PositionType points, VectorPoint& vector_point_field, + VectorGrid vector_grid_field, const VectorHalo vector_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a vector grid value to the points. - Kokkos::deep_copy( vector_point_field, 0.0 ); auto vector_value_g2p = createVectorValueG2P( vector_point_field, -0.5 ); g2p( *vector_grid_field, *vector_halo, points, num_point, Spline<1>(), vector_value_g2p ); - Kokkos::deep_copy( vector_point_host, vector_point_field ); - for ( int p = 0; p < num_point; ++p ) - for ( int d = 0; d < 3; ++d ) +} + +template +void checkVectorValueG2P( const VectorPoint vector_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) + for ( std::size_t d = 0; d < 3; ++d ) EXPECT_FLOAT_EQ( vector_point_host( p, d ), -1.75 ); +} + +template +void scalarGradientG2P( const PositionType points, + VectorPoint& vector_point_field, + ScalarGrid scalar_grid_field, + const ScalarHalo scalar_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *scalar_grid_field, 3.5, Own() ); // Interpolate a scalar grid gradient to the points. - Kokkos::deep_copy( vector_point_field, 0.0 ); auto scalar_gradient_g2p = createScalarGradientG2P( vector_point_field, -0.5 ); g2p( *scalar_grid_field, *scalar_halo, points, num_point, Spline<1>(), scalar_gradient_g2p ); - Kokkos::deep_copy( vector_point_host, vector_point_field ); - for ( int p = 0; p < num_point; ++p ) - for ( int d = 0; d < 3; ++d ) +} + +template +void checkScalarGradientG2P( const VectorPoint vector_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) + for ( std::size_t d = 0; d < 3; ++d ) EXPECT_FLOAT_EQ( vector_point_host( p, d ) + 1.0, 1.0 ); +} + +template +void vectorGradientG2P( const PositionType points, + TensorPoint& tensor_point_field, + VectorGrid vector_grid_field, + const VectorHalo vector_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a vector grid gradient to the points. - Kokkos::deep_copy( tensor_point_field, 0.0 ); auto vector_gradient_g2p = createVectorGradientG2P( tensor_point_field, -0.5 ); g2p( *vector_grid_field, *vector_halo, points, num_point, Spline<1>(), vector_gradient_g2p ); - Kokkos::deep_copy( tensor_point_host, tensor_point_field ); - for ( int p = 0; p < num_point; ++p ) - for ( int i = 0; i < 3; ++i ) - for ( int j = 0; j < 3; ++j ) +} + +template +void checkVectorGradientG2P( const TensorPoint tensor_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) + for ( std::size_t i = 0; i < 3; ++i ) + for ( std::size_t j = 0; j < 3; ++j ) EXPECT_FLOAT_EQ( tensor_point_host( p, i, j ) + 1.0, 1.0 ); +} + +template +void vectorDivergenceG2P( const PositionType points, + ScalarPoint& scalar_point_field, + VectorGrid vector_grid_field, + const VectorHalo vector_halo, + const std::size_t num_point ) +{ + ArrayOp::assign( *vector_grid_field, 3.5, Own() ); // Interpolate a vector grid divergence to the points. - Kokkos::deep_copy( scalar_point_field, 0.0 ); auto vector_div_g2p = createVectorDivergenceG2P( scalar_point_field, -0.5 ); g2p( *vector_grid_field, *vector_halo, points, num_point, Spline<1>(), vector_div_g2p ); - Kokkos::deep_copy( scalar_point_host, scalar_point_field ); - for ( int p = 0; p < num_point; ++p ) +} + +template +void checkVectorDivergenceG2P( const ScalarPoint scalar_point_host, + const std::size_t num_point ) +{ + for ( std::size_t p = 0; p < num_point; ++p ) EXPECT_FLOAT_EQ( scalar_point_host( p ) + 1.0, 1.0 ); } +auto createGrid() +{ + // Create the global mesh. + std::array low_corner = { -1.2, 0.1, 1.1 }; + std::array high_corner = { -0.3, 9.5, 2.3 }; + double cell_size = 0.1; + auto global_mesh = + createUniformGlobalMesh( low_corner, high_corner, cell_size ); + + // Create the global grid. + DimBlockPartitioner<3> partitioner; + std::array is_dim_periodic = { true, true, true }; + auto global_grid = createGlobalGrid( MPI_COMM_WORLD, global_mesh, + is_dim_periodic, partitioner ); + return global_grid; +} + +template +void setPositions( PositionType& points, const LocalGridType local_grid, + const int halo_width ) +{ + auto local_mesh = createLocalMesh( *local_grid ); + auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); + + // Create a point in the center of every cell. + Kokkos::parallel_for( + "fill_points", createExecutionPolicy( cell_space, TEST_EXECSPACE() ), + KOKKOS_LAMBDA( const int i, const int j, const int k ) { + int pi = i - halo_width; + int pj = j - halo_width; + int pk = k - halo_width; + int pid = pi + cell_space.extent( Dim::I ) * + ( pj + cell_space.extent( Dim::J ) * pk ); + int idx[3] = { i, j, k }; + double x[3]; + local_mesh.coordinates( Cell(), idx, x ); + points( pid, Dim::I ) = x[Dim::I]; + points( pid, Dim::J ) = x[Dim::J]; + points( pid, Dim::K ) = x[Dim::K]; + } ); +} +//---------------------------------------------------------------------------// +void interpolationViewTest() +{ + // Create a grid local_grid. + int halo_width = 1; + auto global_grid = createGrid(); + auto local_grid = createLocalGrid( global_grid, halo_width ); + auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); + auto num_point = cell_space.size(); + + // Create particles. + Kokkos::View points( + Kokkos::ViewAllocateWithoutInitializing( "points" ), num_point ); + setPositions( points, local_grid, halo_width ); + + // Node space. + auto node_space = local_grid->indexSpace( Own(), Node(), Local() ); + + // Create a scalar field on the grid. + auto scalar_layout = createArrayLayout( local_grid, 1, Node() ); + auto scalar_grid_field = + createArray( "scalar_grid_field", scalar_layout ); + + // Create a vector field on the grid. + auto vector_layout = createArrayLayout( local_grid, 3, Node() ); + auto vector_grid_field = + createArray( "vector_grid_field", vector_layout ); + + // Create a scalar, vector, and tensor point fields. + Kokkos::View scalar_point_field( + Kokkos::ViewAllocateWithoutInitializing( "scalar_point_field" ), + num_point ); + Kokkos::View vector_point_field( + Kokkos::ViewAllocateWithoutInitializing( "vector_point_field" ), + num_point ); + Kokkos::View tensor_point_field( + Kokkos::ViewAllocateWithoutInitializing( "tensor_point_field" ), + num_point ); + + // Create halos. + auto scalar_halo = + createHalo( NodeHaloPattern<3>(), halo_width, *scalar_grid_field ); + auto vector_halo = + createHalo( NodeHaloPattern<3>(), halo_width, *vector_grid_field ); + + Kokkos::deep_copy( scalar_point_field, 3.5 ); + Kokkos::deep_copy( vector_point_field, 3.5 ); + Kokkos::deep_copy( tensor_point_field, 3.5 ); + + // Check all p2g together. + checkP2G( points, scalar_point_field, vector_point_field, + tensor_point_field, scalar_grid_field, vector_grid_field, + node_space, scalar_halo, vector_halo, num_point ); + + // Separated G2P because of reset/copies (different for AoSoA/View). + Kokkos::deep_copy( scalar_point_field, 0.0 ); + scalarValueG2P( points, scalar_point_field, scalar_grid_field, scalar_halo, + num_point ); + auto scalar_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_point_field ); + checkScalarValueG2P( scalar_point_host, num_point ); + + Kokkos::deep_copy( vector_point_field, 0.0 ); + vectorValueG2P( points, vector_point_field, vector_grid_field, vector_halo, + num_point ); + auto vector_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_point_field ); + checkVectorValueG2P( vector_point_host, num_point ); + + Kokkos::deep_copy( vector_point_field, 0.0 ); + scalarGradientG2P( points, vector_point_field, scalar_grid_field, + scalar_halo, num_point ); + vector_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), vector_point_field ); + checkScalarGradientG2P( vector_point_host, num_point ); + + Kokkos::deep_copy( tensor_point_field, 0.0 ); + vectorGradientG2P( points, tensor_point_field, vector_grid_field, + vector_halo, num_point ); + auto tensor_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), tensor_point_field ); + checkVectorGradientG2P( tensor_point_host, num_point ); + + Kokkos::deep_copy( scalar_point_field, 0.0 ); + vectorDivergenceG2P( points, scalar_point_field, vector_grid_field, + vector_halo, num_point ); + scalar_point_host = Kokkos::create_mirror_view_and_copy( + Kokkos::HostSpace(), scalar_point_field ); + checkVectorDivergenceG2P( scalar_point_host, num_point ); +} + +//---------------------------------------------------------------------------// +void interpolationSliceTest() +{ + // Create a grid local_grid. + int halo_width = 1; + auto global_grid = createGrid(); + auto local_grid = createLocalGrid( global_grid, halo_width ); + auto cell_space = local_grid->indexSpace( Own(), Cell(), Local() ); + auto num_point = cell_space.size(); + + // Create particles. + using member_types = + Cabana::MemberTypes; + Cabana::AoSoA particles( "particles", + num_point ); + auto points = Cabana::slice<2>( particles, "position" ); + setPositions( points, local_grid, halo_width ); + + // Node space. + auto node_space = local_grid->indexSpace( Own(), Node(), Local() ); + + // Create a scalar field on the grid. + auto scalar_layout = createArrayLayout( local_grid, 1, Node() ); + auto scalar_grid_field = + createArray( "scalar_grid_field", scalar_layout ); + + // Create a vector field on the grid. + auto vector_layout = createArrayLayout( local_grid, 3, Node() ); + auto vector_grid_field = + createArray( "vector_grid_field", vector_layout ); + + // Create a scalar, vector, and tensor point fields. + auto tensor_point_field = Cabana::slice<0>( particles, "tensor" ); + auto vector_point_field = Cabana::slice<1>( particles, "vector" ); + auto scalar_point_field = Cabana::slice<3>( particles, "scalar" ); + + // Create halos. + auto scalar_halo = + createHalo( NodeHaloPattern<3>(), halo_width, *scalar_grid_field ); + auto vector_halo = + createHalo( NodeHaloPattern<3>(), halo_width, *vector_grid_field ); + + Cabana::deep_copy( scalar_point_field, 3.5 ); + Cabana::deep_copy( vector_point_field, 3.5 ); + Cabana::deep_copy( tensor_point_field, 3.5 ); + + // Check all p2g together. + checkP2G( points, scalar_point_field, vector_point_field, + tensor_point_field, scalar_grid_field, vector_grid_field, + node_space, scalar_halo, vector_halo, num_point ); + + // Separated to do the reset/copies (different for AoSoA/View). + auto particles_host = + Cabana::create_mirror_view( Kokkos::HostSpace(), particles ); + auto tensor_point_host = Cabana::slice<0>( particles, "tensor" ); + auto vector_point_host = Cabana::slice<1>( particles, "vector" ); + auto scalar_point_host = Cabana::slice<3>( particles, "scalar" ); + + Cabana::deep_copy( scalar_point_field, 0.0 ); + scalarValueG2P( points, scalar_point_field, scalar_grid_field, scalar_halo, + num_point ); + Cabana::deep_copy( scalar_point_host, scalar_point_field ); + checkScalarValueG2P( scalar_point_host, num_point ); + + Cabana::deep_copy( vector_point_field, 0.0 ); + vectorValueG2P( points, vector_point_field, vector_grid_field, vector_halo, + num_point ); + Cabana::deep_copy( vector_point_host, vector_point_field ); + checkVectorValueG2P( vector_point_host, num_point ); + + Cabana::deep_copy( vector_point_field, 0.0 ); + scalarGradientG2P( points, vector_point_field, scalar_grid_field, + scalar_halo, num_point ); + Cabana::deep_copy( vector_point_host, vector_point_field ); + checkScalarGradientG2P( vector_point_host, num_point ); + + Cabana::deep_copy( tensor_point_field, 0.0 ); + vectorGradientG2P( points, tensor_point_field, vector_grid_field, + vector_halo, num_point ); + Cabana::deep_copy( tensor_point_host, tensor_point_field ); + checkVectorGradientG2P( tensor_point_host, num_point ); + + Cabana::deep_copy( scalar_point_field, 0.0 ); + vectorDivergenceG2P( points, scalar_point_field, vector_grid_field, + vector_halo, num_point ); + Cabana::deep_copy( scalar_point_host, scalar_point_field ); + checkVectorDivergenceG2P( scalar_point_host, num_point ); +} + //---------------------------------------------------------------------------// // RUN TESTS //---------------------------------------------------------------------------// -TEST( interpolation, interpolation_test ) { interpolationTest(); } +TEST( interpolation, interpolation_view_test ) { interpolationViewTest(); } + +TEST( interpolation, interpolation_slice_test ) { interpolationSliceTest(); } //---------------------------------------------------------------------------//