From f71924e96fced1e620121f3af4ca65ae7ebd00fa Mon Sep 17 00:00:00 2001 From: Yimin Jin Date: Wed, 12 Feb 2025 16:24:21 +0800 Subject: [PATCH] generalize the portable MF methods to fe_degree < n_q_points_1d cases --- .../matrix_free/portable_evaluation_kernels.h | 620 +++++++++++++ .../matrix_free/portable_fe_evaluation.h | 195 ++-- .../matrix_free/portable_matrix_free.h | 43 +- .../portable_matrix_free.templates.h | 91 +- .../portable_tensor_product_kernels.h | 832 ++++++------------ tests/matrix_free_kokkos/evaluate_1d_shape.cc | 30 +- .../evaluate_1d_shape.output | 24 + tests/matrix_free_kokkos/evaluate_2d_shape.cc | 45 +- .../evaluate_2d_shape.output | 44 +- tests/tests.h | 10 +- 10 files changed, 1219 insertions(+), 715 deletions(-) create mode 100644 include/deal.II/matrix_free/portable_evaluation_kernels.h diff --git a/include/deal.II/matrix_free/portable_evaluation_kernels.h b/include/deal.II/matrix_free/portable_evaluation_kernels.h new file mode 100644 index 0000000000..42a477b42d --- /dev/null +++ b/include/deal.II/matrix_free/portable_evaluation_kernels.h @@ -0,0 +1,620 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2017 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +#ifndef dealii__evaluation_kernels_h +#define dealii__evaluation_kernels_h + +#include + +#include + +#include +#include +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace Portable +{ + namespace internal + { + /** + * Helper function to specify whether a transformation to collocation should + * be used: It should give correct results (first condition), we need to be + * able to initialize the fields in shape_info.templates.h from the + * polynomials (second condition), and it should be the most efficient + * choice in terms of operation counts (third condition). + */ + constexpr bool + use_collocation_evaluation(const unsigned int fe_degree, + const unsigned int n_q_points_1d) + { + // TODO: are the conditions suit for GPU parallelization? + return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) && + (n_q_points_1d <= 3 * fe_degree / 2 + 1); + } + + + + /** + * This struct performs the evaluation of function values and gradients for + * tensor-product finite elements. There are two specialized implementation + * classes FEEvaluationImplCollocation (for Gauss-Lobatto elements where the + * nodal points and the quadrature points coincide and the 'values' + * operation is identity) and FEEvaluationImplTransformToCollocation (which + * can be transformed to a collocation space and can then use the identity + * in these spaces), which both allow for shorter code. + */ + template + struct FEEvaluationImpl + { + using TeamHandle = Kokkos::TeamPolicy< + MemorySpace::Default::kokkos_space::execution_space>::member_type; + using SharedView = Kokkos::View>; + + DEAL_II_HOST_DEVICE static void + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const typename MatrixFree::Data *data, + SharedData *shared_data) + { + if (evaluation_flag == EvaluationFlags::nothing) + return; + + // the evaluator does not need temporary storage since no in-place + // operation takes place in this function + auto scratch_for_eval = + Kokkos::subview(shared_data->scratch_pad, Kokkos::make_pair(0, 0)); + EvaluatorTensorProduct + eval(shared_data->team_member, + data->shape_values, + data->shape_gradients, + data->co_shape_gradients, + scratch_for_eval); + + for (unsigned int c = 0; c < n_components; ++c) + { + auto u = Kokkos::subview(shared_data->values, Kokkos::ALL, c); + auto grad_u = Kokkos::subview(shared_data->gradients, + Kokkos::ALL, + Kokkos::ALL, + c); + + if constexpr (dim == 1) + { + auto temp = + Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, n_q_points_1d)); + + if (evaluation_flag & EvaluationFlags::gradients) + eval.template gradients<0, true, false, false>( + u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); + if (evaluation_flag & EvaluationFlags::values) + { + eval.template values<0, true, false, false>(u, temp); + populate_view(shared_data->team_member, + u, + temp, + n_q_points_1d); + } + } + else if constexpr (dim == 2) + { + constexpr int temp_size = (fe_degree + 1) * n_q_points_1d; + auto temp = Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, temp_size)); + + // grad x + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template gradients<0, true, false, false>(u, temp); + eval.template values<1, true, false, false>( + temp, Kokkos::subview(grad_u, Kokkos::ALL, 0)); + } + + // grad y + eval.template values<0, true, false, false>(u, temp); + if (evaluation_flag & EvaluationFlags::gradients) + eval.template gradients<1, true, false, false>( + temp, Kokkos::subview(grad_u, Kokkos::ALL, 1)); + + // val: can use values applied in x + if (evaluation_flag & EvaluationFlags::values) + eval.template values<1, true, false, false>(temp, u); + } + else if constexpr (dim == 3) + { + constexpr int temp1_size = Utilities::pow(fe_degree + 1, 2) * + n_q_points_1d, + temp2_size = Utilities::pow(n_q_points_1d, 2) * + (fe_degree + 1); + + auto temp1 = Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, temp1_size)); + auto temp2 = + Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(temp1_size, + temp1_size + temp2_size)); + + if (evaluation_flag & EvaluationFlags::gradients) + { + // grad x + eval.template gradients<0, true, false, false>(u, temp1); + eval.template values<1, true, false, false>(temp1, temp2); + eval.template values<2, true, false, false>( + temp2, Kokkos::subview(grad_u, Kokkos::ALL, 0)); + } + + // grad y + eval.template values<0, true, false, false>(u, temp1); + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template gradients<1, true, false, false>(temp1, + temp2); + eval.template values<2, true, false, false>( + temp2, Kokkos::subview(grad_u, Kokkos::ALL, 1)); + } + + // grad z: can use the values applied in x direction stored + // in temp1 + eval.template values<1, true, false, false>(temp1, temp2); + if (evaluation_flag & EvaluationFlags::gradients) + eval.template gradients<2, true, false, false>( + temp2, Kokkos::subview(grad_u, Kokkos::ALL, 2)); + + // val: can use the values applied in x & y direction stored + // in temp2 + if (evaluation_flag & EvaluationFlags::values) + eval.template values<2, true, false, false>(temp2, u); + } + else + Assert(false, ExcMessage("dim must not exceed 3!")); + } + } + + + + DEAL_II_HOST_DEVICE static void + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const typename MatrixFree::Data *data, + SharedData *shared_data) + { + if (integration_flag == EvaluationFlags::nothing) + return; + + // the evaluator does not need temporary storage since no in-place + // operation takes place in this function + auto scratch_for_eval = + Kokkos::subview(shared_data->scratch_pad, Kokkos::make_pair(0, 0)); + EvaluatorTensorProduct + eval(shared_data->team_member, + data->shape_values, + data->shape_gradients, + data->co_shape_gradients, + scratch_for_eval); + + for (unsigned int c = 0; c < n_components; ++c) + { + auto u = Kokkos::subview(shared_data->values, Kokkos::ALL, c); + auto grad_u = Kokkos::subview(shared_data->gradients, + Kokkos::ALL, + Kokkos::ALL, + c); + + if constexpr (dim == 1) + { + auto temp = + Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, fe_degree + 1)); + + if ((integration_flag & EvaluationFlags::values) && + !(integration_flag & EvaluationFlags::gradients)) + { + eval.template values<0, false, false, false>(u, temp); + populate_view(shared_data->team_member, + u, + temp, + fe_degree + 1); + } + if (integration_flag & EvaluationFlags::gradients) + { + if (integration_flag & EvaluationFlags::values) + { + eval.template values<0, false, false, false>(u, temp); + eval.template gradients<0, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), temp); + populate_view(shared_data->team_member, + u, + temp, + fe_degree + 1); + } + else + eval.template gradients<0, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), u); + } + } + else if constexpr (dim == 2) + { + constexpr int temp_size = (fe_degree + 1) * n_q_points_1d; + auto temp = Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, temp_size)); + + if ((integration_flag & EvaluationFlags::values) && + !(integration_flag & EvaluationFlags::gradients)) + { + eval.template values<1, false, false, false>(u, temp); + eval.template values<0, false, false, false>(temp, u); + } + if (integration_flag & EvaluationFlags::gradients) + { + eval.template gradients<1, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), temp); + if (integration_flag & EvaluationFlags::values) + eval.template values<1, false, true, false>(u, temp); + eval.template values<0, false, false, false>(temp, u); + eval.template values<1, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), temp); + eval.template gradients<0, false, true, false>(temp, u); + } + } + else if constexpr (dim == 3) + { + constexpr int temp1_size = Utilities::pow(n_q_points_1d, 2) * + (fe_degree + 1), + temp2_size = Utilities::pow(fe_degree + 1, 2) * + n_q_points_1d; + + auto temp1 = Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, temp1_size)); + auto temp2 = + Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(temp1_size, + temp1_size + temp2_size)); + + if ((integration_flag & EvaluationFlags::values) && + !(integration_flag & EvaluationFlags::gradients)) + { + eval.template values<2, false, false, false>(u, temp1); + eval.template values<1, false, false, false>(temp1, temp2); + eval.template values<0, false, false, false>(temp2, u); + } + if (integration_flag & EvaluationFlags::gradients) + { + eval.template gradients<2, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 2), temp1); + if (integration_flag & EvaluationFlags::values) + eval.template values<2, false, true, false>(u, temp1); + eval.template values<1, false, false, false>(temp1, temp2); + eval.template values<2, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), temp1); + eval.template gradients<1, false, true, false>(temp1, + temp2); + eval.template values<0, false, false, false>(temp2, u); + eval.template values<2, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), temp1); + eval.template values<1, false, false, false>(temp1, temp2); + eval.template gradients<0, false, true, false>(temp2, u); + } + } + else + Assert(false, ExcMessage("dim must not exceed 3!")); + } + } + }; + + + + /** + * This struct performs the evaluation of function values and gradients for + * tensor-product finite elements. This is a specialization for elements + * where the nodal points coincide with the quadrature points like FE_Q + * shape functions on Gauss-Lobatto elements integrated with Gauss-Lobatto + * quadrature. The assumption of this class is that the shape 'values' + * operation is identity, which allows us to write shorter code. + * + * In literature, this form of evaluation is often called spectral + * evaluation, spectral collocation or simply collocation, meaning the same + * location for shape functions and evaluation space (quadrature points). + */ + template + struct FEEvaluationImplCollocation + { + DEAL_II_HOST_DEVICE static void + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const typename MatrixFree::Data *data, + SharedData *shared_data) + { + // since the dof values have already been stored in + // shared_data->values, there is nothing to do if the gradients are + // not required + if (!(evaluation_flag & EvaluationFlags::gradients)) + return; + + constexpr int n_points = Utilities::pow(fe_degree + 1, dim); + auto scratch_for_eval = Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, n_points)); + + EvaluatorTensorProduct + eval(shared_data->team_member, + data->shape_values, + data->shape_gradients, + data->co_shape_gradients, + scratch_for_eval); + + for (unsigned int c = 0; c < n_components; ++c) + { + auto u = Kokkos::subview(shared_data->values, Kokkos::ALL, c); + auto grad_u = Kokkos::subview(shared_data->gradients, + Kokkos::ALL, + Kokkos::ALL, + c); + + eval.template co_gradients<0, true, false, false>( + u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); + if constexpr (dim > 1) + eval.template co_gradients<1, true, false, false>( + u, Kokkos::subview(grad_u, Kokkos::ALL, 1)); + if constexpr (dim > 2) + eval.template co_gradients<2, true, false, false>( + u, Kokkos::subview(grad_u, Kokkos::ALL, 2)); + } + } + + + DEAL_II_HOST_DEVICE static void + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const typename MatrixFree::Data *data, + SharedData *shared_data) + { + // since the quad values have already been stored in + // shared_data->values, there is nothing to do if the gradients are + // not required + if (!(integration_flag & EvaluationFlags::gradients)) + return; + + constexpr int n_points = Utilities::pow(fe_degree + 1, dim); + auto scratch_for_eval = Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, n_points)); + + EvaluatorTensorProduct + eval(shared_data->team_member, + data->shape_values, + data->shape_gradients, + data->co_shape_gradients, + scratch_for_eval); + + for (unsigned int c = 0; c < n_components; ++c) + { + auto u = Kokkos::subview(shared_data->values, Kokkos::ALL, c); + auto grad_u = Kokkos::subview(shared_data->gradients, + Kokkos::ALL, + Kokkos::ALL, + c); + + if constexpr (dim == 1) + { + if (integration_flag & EvaluationFlags::values) + eval.template co_gradients<0, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), u); + else + eval.template co_gradients<2, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 2), u); + } + else if constexpr (dim == 2) + { + if (integration_flag & EvaluationFlags::values) + eval.template co_gradients<1, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), u); + else + eval.template co_gradients<1, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), u); + eval.template co_gradients<0, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), u); + } + else if constexpr (dim == 3) + { + if (integration_flag & EvaluationFlags::values) + eval.template co_gradients<2, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 2), u); + else + eval.template co_gradients<2, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 2), u); + eval.template co_gradients<1, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), u); + eval.template co_gradients<0, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), u); + } + else + Assert(false, ExcMessage("dim must not exceed 3!")); + } + } + }; + + + + /** + * This struct performs the evaluation of function values and gradients for + * tensor-product finite elements. This is a specialization for symmetric + * basis functions about the mid point 0.5 of the unit interval with the + * same number of quadrature points as degrees of freedom. In that case, we + * can first transform the basis to one that has the nodal points in the + * quadrature points (i.e., the collocation space) and then perform the + * evaluation of the first and second derivatives in this transformed space, + * using the identity operation for the shape values. + */ + template + struct FEEvaluationImplTransformToCollocation + { + DEAL_II_HOST_DEVICE static void + evaluate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const typename MatrixFree::Data *data, + SharedData *shared_data) + { + constexpr int scratch_size = Utilities::pow(n_q_points_1d, dim); + auto scratch_for_eval = + Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, scratch_size)); + + EvaluatorTensorProduct + eval(shared_data->team_member, + data->shape_values, + data->shape_gradients, + data->co_shape_gradients, + scratch_for_eval); + + for (unsigned int c = 0; c < n_components; ++c) + { + auto u = Kokkos::subview(shared_data->values, Kokkos::ALL, c); + auto grad_u = Kokkos::subview(shared_data->gradients, + Kokkos::ALL, + Kokkos::ALL, + c); + + eval.template values<0, true, false, true>(u, u); + if constexpr (dim > 1) + eval.template values<1, true, false, true>(u, u); + if constexpr (dim > 2) + eval.template values<2, true, false, true>(u, u); + + if (evaluation_flag & EvaluationFlags::gradients) + { + eval.template co_gradients<0, true, false, false>( + u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); + if constexpr (dim > 1) + eval.template co_gradients<1, true, false, false>( + u, Kokkos::subview(grad_u, Kokkos::ALL, 1)); + if constexpr (dim > 2) + eval.template co_gradients<2, true, false, false>( + u, Kokkos::subview(grad_u, Kokkos::ALL, 2)); + } + } + } + + + DEAL_II_HOST_DEVICE static void + integrate(const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + const typename MatrixFree::Data *data, + const SharedData *shared_data) + { + constexpr int scratch_size = Utilities::pow(n_q_points_1d, dim); + auto scratch_for_eval = + Kokkos::subview(shared_data->scratch_pad, + Kokkos::make_pair(0, scratch_size)); + + EvaluatorTensorProduct + eval(shared_data->team_member, + data->shape_values, + data->shape_gradients, + data->co_shape_gradients, + scratch_for_eval); + + for (unsigned int c = 0; c < n_components; ++c) + { + auto u = Kokkos::subview(shared_data->values, Kokkos::ALL, c); + auto grad_u = Kokkos::subview(shared_data->gradients, + Kokkos::ALL, + Kokkos::ALL, + c); + + // apply derivatives in collocation space + if (integration_flag & EvaluationFlags::gradients) + { + if constexpr (dim == 1) + { + if (integration_flag & EvaluationFlags::values) + eval.template co_gradients<0, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), u); + else + eval.template co_gradients<2, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 2), u); + } + else if constexpr (dim == 2) + { + if (integration_flag & EvaluationFlags::values) + eval.template co_gradients<1, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), u); + else + eval.template co_gradients<1, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), u); + eval.template co_gradients<0, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), u); + } + else if constexpr (dim == 3) + { + if (integration_flag & EvaluationFlags::values) + eval.template co_gradients<2, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 2), u); + else + eval.template co_gradients<2, false, false, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 2), u); + eval.template co_gradients<1, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 1), u); + eval.template co_gradients<0, false, true, false>( + Kokkos::subview(grad_u, Kokkos::ALL, 0), u); + } + else + Assert(false, ExcMessage("dim must not exceed 3!")); + } + + // transform back to the original space + if constexpr (dim > 2) + eval.template values<2, false, false, true>(u, u); + if constexpr (dim > 1) + eval.template values<1, false, false, true>(u, u); + eval.template values<0, false, false, true>(u, u); + } + } + }; + } // end of namespace internal +} // end of namespace Portable + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/portable_fe_evaluation.h b/include/deal.II/matrix_free/portable_fe_evaluation.h index dfa5a46ff6..ae23f70eb8 100644 --- a/include/deal.II/matrix_free/portable_fe_evaluation.h +++ b/include/deal.II/matrix_free/portable_fe_evaluation.h @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -106,10 +107,18 @@ namespace Portable Utilities::pow(n_q_points_1d, dim); /** - * Number of tensor degrees of freedoms per cell. + * Number of tensor degrees of freedom of a scalar component determined + * from the given template argument `fe_degree`. + */ + static constexpr unsigned int tensor_dofs_per_component = + Utilities::pow(fe_degree + 1, dim); + + /** + * Number of tensor degrees of freedom of all component determined from + * the given template argument `fe_degree`. */ static constexpr unsigned int tensor_dofs_per_cell = - Utilities::pow(fe_degree + 1, dim) * n_components; + tensor_dofs_per_component * n_components; /** * Constructor. @@ -271,13 +280,14 @@ namespace Portable read_dof_values(const Number *src) { // Populate the scratch memory - Kokkos::parallel_for( - Kokkos::TeamThreadRange(shared_data->team_member, n_q_points), - [&](const int &i) { - for (unsigned int c = 0; c < n_components_; ++c) - shared_data->values(i, c) = src[data->local_to_global( - cell_id, i + (tensor_dofs_per_cell / n_components) * c)]; - }); + Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member, + tensor_dofs_per_component), + [&](const int &i) { + for (unsigned int c = 0; c < n_components_; ++c) + shared_data->values(i, c) = + src[data->local_to_global( + cell_id, i + tensor_dofs_per_component * c)]; + }); shared_data->team_member.team_barrier(); for (unsigned int c = 0; c < n_components_; ++c) @@ -313,24 +323,25 @@ namespace Portable if (data->use_coloring) { Kokkos::parallel_for( - Kokkos::TeamThreadRange(shared_data->team_member, n_q_points), + Kokkos::TeamThreadRange(shared_data->team_member, + tensor_dofs_per_component), [&](const int &i) { for (unsigned int c = 0; c < n_components_; ++c) - dst[data->local_to_global( - cell_id, i + (tensor_dofs_per_cell / n_components) * c)] += + dst[data->local_to_global(cell_id, + i + tensor_dofs_per_component * c)] += shared_data->values(i, c); }); } else { Kokkos::parallel_for( - Kokkos::TeamThreadRange(shared_data->team_member, n_q_points), + Kokkos::TeamThreadRange(shared_data->team_member, + tensor_dofs_per_component), [&](const int &i) { for (unsigned int c = 0; c < n_components_; ++c) - Kokkos::atomic_add( - &dst[data->local_to_global( - cell_id, i + (tensor_dofs_per_cell / n_components) * c)], - shared_data->values(i, c)); + Kokkos::atomic_add(&dst[data->local_to_global( + cell_id, i + (tensor_dofs_per_component)*c)], + shared_data->values(i, c)); }); } } @@ -340,50 +351,42 @@ namespace Portable template DEAL_II_HOST_DEVICE void - FEEvaluation::evaluate( - const EvaluationFlags::EvaluationFlags evaluate_flag) + FEEvaluation::evaluate( + const EvaluationFlags::EvaluationFlags evaluation_flag) { - // First evaluate the gradients because it requires values that will be - // changed if evaluate_val is true - internal::EvaluatorTensorProduct< - internal::EvaluatorVariant::evaluate_general, - dim, - fe_degree, - n_q_points_1d, - Number> - evaluator_tensor_product(shared_data->team_member, - data->shape_values, - data->shape_gradients, - data->co_shape_gradients); + using ElementType = ::dealii::internal::MatrixFreeFunctions::ElementType; - for (unsigned int c = 0; c < n_components_; ++c) + if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && + data->element_type == ElementType::tensor_symmetric_collocation) { - if ((evaluate_flag & EvaluationFlags::values) && - (evaluate_flag & EvaluationFlags::gradients)) - { - evaluator_tensor_product.evaluate_values_and_gradients( - Kokkos::subview(shared_data->values, Kokkos::ALL, c), - Kokkos::subview( - shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c)); - shared_data->team_member.team_barrier(); - } - else if (evaluate_flag & EvaluationFlags::gradients) - { - evaluator_tensor_product.evaluate_gradients( - Kokkos::subview(shared_data->values, Kokkos::ALL, c), - Kokkos::subview( - shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c)); - shared_data->team_member.team_barrier(); - } - else if (evaluate_flag & EvaluationFlags::values) - { - evaluator_tensor_product.evaluate_values( - Kokkos::subview(shared_data->values, Kokkos::ALL, c)); - shared_data->team_member.team_barrier(); - } + internal::FEEvaluationImplCollocation::evaluate( + n_components, evaluation_flag, data, shared_data); + } + // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see + // shape_info.h for more details + else if (fe_degree >= 0 && + internal::use_collocation_evaluation(fe_degree, n_q_points_1d) && + data->element_type <= ElementType::tensor_symmetric) + { + internal::FEEvaluationImplTransformToCollocation< + dim, + fe_degree, + n_q_points_1d, + Number>::evaluate(n_components, evaluation_flag, data, shared_data); + } + else if (fe_degree >= 0 && + data->element_type <= ElementType::tensor_symmetric_no_collocation) + { + internal::FEEvaluationImpl:: + evaluate(n_components, evaluation_flag, data, shared_data); + } + else + { + Kokkos::abort("The element type is not yet supported by the portable " + "matrix-free module."); } } @@ -415,42 +418,36 @@ namespace Portable FEEvaluation::integrate( const EvaluationFlags::EvaluationFlags integration_flag) { - internal::EvaluatorTensorProduct< - internal::EvaluatorVariant::evaluate_general, - dim, - fe_degree, - n_q_points_1d, - Number> - evaluator_tensor_product(shared_data->team_member, - data->shape_values, - data->shape_gradients, - data->co_shape_gradients); - + using ElementType = ::dealii::internal::MatrixFreeFunctions::ElementType; - for (unsigned int c = 0; c < n_components_; ++c) + if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && + data->element_type == ElementType::tensor_symmetric_collocation) { - if ((integration_flag & EvaluationFlags::values) && - (integration_flag & EvaluationFlags::gradients)) - { - evaluator_tensor_product.integrate_values_and_gradients( - Kokkos::subview(shared_data->values, Kokkos::ALL, c), - Kokkos::subview( - shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c)); - } - else if (integration_flag & EvaluationFlags::values) - { - evaluator_tensor_product.integrate_values( - Kokkos::subview(shared_data->values, Kokkos::ALL, c)); - shared_data->team_member.team_barrier(); - } - else if (integration_flag & EvaluationFlags::gradients) - { - evaluator_tensor_product.template integrate_gradients( - Kokkos::subview(shared_data->values, Kokkos::ALL, c), - Kokkos::subview( - shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c)); - shared_data->team_member.team_barrier(); - } + internal::FEEvaluationImplCollocation:: + integrate(n_components, integration_flag, data, shared_data); + } + // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see + // shape_info.h for more details + else if (fe_degree >= 0 && + internal::use_collocation_evaluation(fe_degree, n_q_points_1d) && + data->element_type <= ElementType::tensor_symmetric) + { + internal::FEEvaluationImplTransformToCollocation< + dim, + fe_degree, + n_q_points_1d, + Number>::integrate(n_components, integration_flag, data, shared_data); + } + else if (fe_degree >= 0 && + data->element_type <= ElementType::tensor_symmetric_no_collocation) + { + internal::FEEvaluationImpl:: + integrate(n_components, integration_flag, data, shared_data); + } + else + { + Kokkos::abort("The element type is not yet supported by the portable " + "matrix-free module."); } } @@ -486,6 +483,7 @@ namespace Portable FEEvaluation::get_value( int q_point) const { + Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError()); if constexpr (n_components_ == 1) { return shared_data->values(q_point, 0); @@ -512,17 +510,18 @@ namespace Portable n_components_, Number>::value_type FEEvaluation:: - get_dof_value(int q_point) const + get_dof_value(int dof) const { + Assert(dof >= 0 && dof < tensor_dofs_per_component, ExcInternalError()); if constexpr (n_components_ == 1) { - return shared_data->values(q_point, 0); + return shared_data->values(dof, 0); } else { value_type result; for (unsigned int c = 0; c < n_components; ++c) - result[c] = shared_data->values(q_point, c); + result[c] = shared_data->values(dof, c); return result; } } @@ -538,6 +537,7 @@ namespace Portable FEEvaluation:: submit_value(const value_type &val_in, int q_point) { + Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError()); if constexpr (n_components_ == 1) { shared_data->values(q_point, 0) = val_in * data->JxW(cell_id, q_point); @@ -559,16 +559,17 @@ namespace Portable typename Number> DEAL_II_HOST_DEVICE void FEEvaluation:: - submit_dof_value(const value_type &val_in, int q_point) + submit_dof_value(const value_type &val_in, int dof) { + Assert(dof >= 0 && dof < tensor_dofs_per_component, ExcInternalError()); if constexpr (n_components_ == 1) { - shared_data->values(q_point, 0) = val_in; + shared_data->values(dof, 0) = val_in; } else { for (unsigned int c = 0; c < n_components; ++c) - shared_data->values(q_point, c) = val_in[c]; + shared_data->values(dof, c) = val_in[c]; } } @@ -587,6 +588,7 @@ namespace Portable FEEvaluation:: get_gradient(int q_point) const { + Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError()); gradient_type grad; if constexpr (n_components_ == 1) @@ -627,6 +629,7 @@ namespace Portable FEEvaluation:: submit_gradient(const gradient_type &grad_in, int q_point) { + Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError()); if constexpr (n_components_ == 1) { for (unsigned int d_1 = 0; d_1 < dim; ++d_1) diff --git a/include/deal.II/matrix_free/portable_matrix_free.h b/include/deal.II/matrix_free/portable_matrix_free.h index dc1b9c22cf..dee0088832 100644 --- a/include/deal.II/matrix_free/portable_matrix_free.h +++ b/include/deal.II/matrix_free/portable_matrix_free.h @@ -35,6 +35,8 @@ #include #include +#include + #include @@ -234,6 +236,18 @@ namespace Portable */ bool use_coloring; + /** + * Encodes the type of element detected at construction. FEEvaluation + * will select the most efficient algorithm based on the given element + * type. + */ + ::dealii::internal::MatrixFreeFunctions::ElementType element_type; + + /** + * Size of the scratch pad for temporary storage in shared memory. + */ + unsigned int scratch_pad_size; + /** * Return the quadrature point index local. The index is * only unique for a given MPI process. @@ -474,6 +488,18 @@ namespace Portable */ types::global_dof_index n_dofs; + /** + * Encodes the type of element detected at construction. FEEvaluation + * will select the most efficient algorithm based on the given element + * type. + */ + ::dealii::internal::MatrixFreeFunctions::ElementType element_type; + + /** + * Size of the scratch pad for temporary storage in shared memory. + */ + unsigned int scratch_pad_size; + /** * Degree of the finite element used. */ @@ -627,14 +653,20 @@ namespace Portable Number ***, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits>; + using SharedViewScratchPad = Kokkos::View< + Number *, + MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, + Kokkos::MemoryTraits>; DEAL_II_HOST_DEVICE - SharedData(const TeamHandle &team_member, - const SharedViewValues &values, - const SharedViewGradients &gradients) + SharedData(const TeamHandle &team_member, + const SharedViewValues &values, + const SharedViewGradients &gradients, + const SharedViewScratchPad &scratch_pad) : team_member(team_member) , values(values) , gradients(gradients) + , scratch_pad(scratch_pad) {} /** @@ -651,6 +683,11 @@ namespace Portable * Memory for computed gradients in reference coordinate system. */ SharedViewGradients gradients; + + /** + * Memory for temporary arrays required by evaluation and integration. + */ + SharedViewScratchPad scratch_pad; }; diff --git a/include/deal.II/matrix_free/portable_matrix_free.templates.h b/include/deal.II/matrix_free/portable_matrix_free.templates.h index e8c760923b..3903861aa7 100644 --- a/include/deal.II/matrix_free/portable_matrix_free.templates.h +++ b/include/deal.II/matrix_free/portable_matrix_free.templates.h @@ -345,6 +345,11 @@ namespace Portable MemorySpace::Default::kokkos_space::execution_space:: scratch_memory_space, Kokkos::MemoryTraits>; + using SharedViewScratchPad = + Kokkos::View>; ApplyKernel(Functor func, const typename MatrixFree::Data gpu_data, @@ -367,13 +372,12 @@ namespace Portable size_t team_shmem_size(int /*team_size*/) const { - return SharedViewValues::shmem_size(Functor::n_local_dofs / - gpu_data.n_components, + return SharedViewValues::shmem_size(Functor::n_q_points, gpu_data.n_components) + - SharedViewGradients::shmem_size(Functor::n_local_dofs / - gpu_data.n_components, + SharedViewGradients::shmem_size(Functor::n_q_points, dim, - gpu_data.n_components); + gpu_data.n_components) + + SharedViewScratchPad::shmem_size(gpu_data.scratch_pad_size); } @@ -382,16 +386,20 @@ namespace Portable operator()(const TeamHandle &team_member) const { // Get the scratch memory - SharedViewValues values(team_member.team_shmem(), - Functor::n_local_dofs / gpu_data.n_components, + SharedViewValues values(team_member.team_shmem(), + Functor::n_q_points, gpu_data.n_components); - SharedViewGradients gradients(team_member.team_shmem(), - Functor::n_local_dofs / - gpu_data.n_components, + SharedViewGradients gradients(team_member.team_shmem(), + Functor::n_q_points, dim, gpu_data.n_components); + SharedViewScratchPad scratch_pad(team_member.team_shmem(), + gpu_data.scratch_pad_size); - SharedData shared_data(team_member, values, gradients); + SharedData shared_data(team_member, + values, + gradients, + scratch_pad); func(team_member.league_rank(), &gpu_data, &shared_data, src, dst); } }; @@ -500,6 +508,8 @@ namespace Portable data_copy.padding_length = padding_length; data_copy.row_start = row_start[color]; data_copy.use_coloring = use_coloring; + data_copy.element_type = element_type; + data_copy.scratch_pad_size = scratch_pad_size; return data_copy; } @@ -660,6 +670,49 @@ namespace Portable + namespace + { + /** + * Helper function for determining the scratch pad size. + */ + inline unsigned int + compute_scratch_pad_size( + const ::dealii::internal::MatrixFreeFunctions::ElementType element_type, + const int dim, + const int fe_degree, + const int n_q_points_1d) + { + using ElementType = ::dealii::internal::MatrixFreeFunctions::ElementType; + + if (fe_degree >= 0 && element_type <= ElementType::tensor_symmetric) + { + // evaluate/integrate with FEEvaluationImplCollocation or + // FEEvaluationImplTransformToCollocation + return Utilities::pow(n_q_points_1d, dim); + } + else if (fe_degree >= 0 && + element_type <= ElementType::tensor_symmetric_no_collocation) + { + // evaluate/integrate with FEEvaluationImpl + if (dim == 1) + return n_q_points_1d; + else if (dim == 2) + return (fe_degree + 1) * n_q_points_1d; + else if (dim == 3) + return (fe_degree + 1) * n_q_points_1d * + (fe_degree + 1 + n_q_points_1d); + else + AssertThrow(false, ExcNotImplemented()); + } + else + AssertThrow(false, ExcNotImplemented()); + + return numbers::invalid_unsigned_int; + } + } // namespace + + + template template void @@ -691,13 +744,15 @@ namespace Portable const unsigned int n_dofs_1d = fe_degree + 1; const unsigned int n_q_points_1d = quad.size(); - Assert(n_dofs_1d == n_q_points_1d, - ExcMessage("n_q_points_1d must be equal to fe_degree+1.")); + // TODO remove the limitation in the future + AssertThrow(n_dofs_1d <= n_q_points_1d, + ExcMessage("n_q_points_1d must be greater than or equal to " + "fe_degree + 1.")); // Set padding length to the closest power of two larger than or equal to // the number of threads. - padding_length = 1 << static_cast( - std::ceil(dim * std::log2(fe_degree + 1.))); + padding_length = 1 << static_cast(std::ceil( + dim * std::log2(static_cast(n_q_points_1d)))); dofs_per_cell = fe.n_dofs_per_cell(); n_components = fe.n_components(); @@ -707,6 +762,12 @@ namespace Portable ::dealii::internal::MatrixFreeFunctions::ShapeInfo shape_info(quad, fe); + this->element_type = shape_info.element_type; + this->scratch_pad_size = compute_scratch_pad_size(shape_info.element_type, + dim, + fe_degree, + n_q_points_1d); + unsigned int size_shape_values = n_dofs_1d * n_q_points_1d; shape_values = Kokkos::View( diff --git a/include/deal.II/matrix_free/portable_tensor_product_kernels.h b/include/deal.II/matrix_free/portable_tensor_product_kernels.h index 6e29a3a593..e168512b18 100644 --- a/include/deal.II/matrix_free/portable_tensor_product_kernels.h +++ b/include/deal.II/matrix_free/portable_tensor_product_kernels.h @@ -42,16 +42,46 @@ namespace Portable + /** + * Helper function that copies or adds the first N entries of src to + * dst, depending on the template argument "add". + */ + template + DEAL_II_HOST_DEVICE void + populate_view( + const Kokkos::TeamPolicy< + MemorySpace::Default::kokkos_space::execution_space>::member_type + &team_member, + ViewTypeOut dst, + const ViewTypeIn src, + const int N) + { + Assert(dst.size() >= N, ExcInternalError()); + Assert(src.size() >= N, ExcInternalError()); + + Kokkos::parallel_for(Kokkos::TeamVectorRange(team_member, N), + [&](const int i) { + if constexpr (add) + Kokkos::atomic_add(&dst(i), src(i)); + else + dst(i) = src(i); + }); + + team_member.team_barrier(); + } + + + #if KOKKOS_VERSION >= 40000 /** - * Helper function for values() and gradients() in 1D + * Helper function for apply() in 1D */ - template DEAL_II_HOST_DEVICE void @@ -63,49 +93,44 @@ namespace Portable const ViewTypeIn in, ViewTypeOut out) { - Number t[n_q_points_1d]; - Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, n_q_points_1d), - [&](const int &q) { - t[q] = 0; - // This loop simply multiplies the shape function - // at the quadrature point by the value finite - // element coefficient. - // FIXME check why using parallel_reduce - // ThreadVector is slower - for (int k = 0; k < n_q_points_1d; ++k) + constexpr int Nk = (contract_over_rows ? n_rows : n_columns), + Nq = (contract_over_rows ? n_columns : n_rows); + + Assert(shape_data.size() == n_rows * n_columns, ExcInternalError()); + Assert(in.size() >= Nk, ExcInternalError()); + Assert(out.size() >= Nq, ExcInternalError()); + + Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, Nq), + [&](const int q) { + Number sum = 0; + for (int k = 0; k < Nk; ++k) { - const unsigned int shape_idx = - dof_to_quad ? (q + k * n_q_points_1d) : - (k + q * n_q_points_1d); - const unsigned int source_idx = k; - t[q] += shape_data[shape_idx] * in(source_idx); + const int shape_idx = + (contract_over_rows ? q + k * Nq : + k + q * Nk); + sum += shape_data(shape_idx) * in(k); } - }); - - if constexpr (in_place) - team_member.team_barrier(); - Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, n_q_points_1d), - [&](const int &q) { - const unsigned int destination_idx = q; if constexpr (add) - Kokkos::atomic_add(&out(destination_idx), t[q]); + Kokkos::atomic_add(&out(q), sum); else - out(destination_idx) = t[q]; + out(q) = sum; }); + + team_member.team_barrier(); } /** - * Helper function for values() and gradients() in 2D + * Helper function for apply() in 2D */ - template DEAL_II_HOST_DEVICE void @@ -119,58 +144,70 @@ namespace Portable { using TeamType = Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space>::member_type; - constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 2); - Number t[n_q_points]; - auto thread_policy = + // Sizes of the input and output vectors: + // ----------------------------------------------------------- + // direction | contract_over_rows | !contract_over_rows + // ----------------------------------------------------------- + // 0 | m x m -> n x m | n x m -> m x m + // ----------------------------------------------------------- + // 1 | n x m -> n x n | n x n -> n x m + // ----------------------------------------------------------- + // + // Directions of the cycle indices: + // ----------------------------- + // direction | j | q/k + // ----------------------------- + // 0 | 1 | 0 + // ----------------------------- + // 1 | 0 | 1 + // ----------------------------- + constexpr int Nj = (direction < 1 ? n_rows : n_columns), + Nk = (contract_over_rows ? n_rows : n_columns), + Nq = (contract_over_rows ? n_columns : n_rows); + + Assert(shape_data.size() == n_rows * n_columns, ExcInternalError()); + Assert(in.size() >= Nj * Nk, ExcInternalError()); + Assert(out.size() >= Nj * Nq, ExcInternalError()); + + auto thread_policy = Kokkos::TeamThreadMDRange, TeamType>(team_member, - n_q_points_1d, - n_q_points_1d); - Kokkos::parallel_for(thread_policy, [&](const int i, const int j) { - int q_point = i + j * n_q_points_1d; - - // This loop simply multiplies the shape function at the quadrature - // point by the value finite element coefficient. - // FIXME check why using parallel_reduce ThreadVector is slower - const int base_shape = dof_to_quad ? j : j * n_q_points_1d; - const int stride_shape = dof_to_quad ? n_q_points_1d : 1; - const int base_in = (direction == 0) ? (n_q_points_1d * i) : i; - const int stride_in = Utilities::pow(n_q_points_1d, direction); - Number sum = shape_data[base_shape] * in(base_in); - for (int k = 1; k < n_q_points_1d; ++k) - { - sum += shape_data[base_shape + k * stride_shape] * - in(base_in + k * stride_in); - } - t[q_point] = sum; - }); + Nj, + Nq); + Kokkos::parallel_for(thread_policy, [&](const int j, const int q) { + const int base_shape = contract_over_rows ? q : q * n_columns; + const int stride_shape = contract_over_rows ? n_columns : 1; - if (in_place) - team_member.team_barrier(); + const int base_in = (direction == 0 ? j * Nk : j); + const int stride_in = Utilities::pow(n_columns, direction); - Kokkos::parallel_for(thread_policy, [&](const int i, const int j) { - const int q_point = i + j * n_q_points_1d; - const unsigned int destination_idx = - (direction == 0) ? (j + n_q_points_1d * i) : (i + n_q_points_1d * j); + Number sum = shape_data(base_shape) * in(base_in); + for (int k = 1; k < Nk; ++k) + sum += shape_data(base_shape + k * stride_shape) * + in(base_in + k * stride_in); - if (add) - Kokkos::atomic_add(&out(destination_idx), t[q_point]); + const int index_out = (direction == 0 ? j * Nq + q : j + q * Nj); + + if constexpr (add) + Kokkos::atomic_add(&out(index_out), sum); else - out(destination_idx) = t[q_point]; + out(index_out) = sum; }); + + team_member.team_barrier(); } /** - * Helper function for values() and gradients() in 3D + * Helper function for apply() in 3D */ - template DEAL_II_HOST_DEVICE void @@ -184,68 +221,78 @@ namespace Portable { using TeamType = Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space>::member_type; - constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 3); - Number t[n_q_points]; + // Sizes of the input and output vectors: + // ------------------------------------------------------------------ + // direction | contract_over_rows | !contract_over_rows + // ------------------------------------------------------------------ + // 0 | m x m x m -> n x m x m | n x m x m -> m x m x m + // ------------------------------------------------------------------ + // 1 | n x m x m -> n x n x m | n x n x m -> n x m x m + // ------------------------------------------------------------------ + // 2 | n x n x m -> n x n x n | n x n x n -> n x n x m + // ------------------------------------------------------------------ + // + // Directions of the cycle indices: + // ------------------------------------- + // direction | i | j | q/k + // ------------------------------------- + // 0 | 2 | 1 | 0 + // ------------------------------------- + // 1 | 0 | 2 | 1 + // ------------------------------------- + // 2 | 1 | 0 | 2 + // ------------------------------------- + constexpr int Ni = (direction < 1 ? n_rows : n_columns), + Nj = (direction < 2 ? n_rows : n_columns), + Nk = (contract_over_rows ? n_rows : n_columns), + Nq = (contract_over_rows ? n_columns : n_rows); + + Assert(shape_data.size() == n_rows * n_columns, ExcInternalError()); + Assert(in.size() >= Ni * Nj * Nk, ExcInternalError()); + Assert(out.size() >= Ni * Nj * Nq, ExcInternalError()); + auto thread_policy = Kokkos::TeamThreadMDRange, TeamType>( - team_member, n_q_points_1d, n_q_points_1d, n_q_points_1d); + team_member, Ni, Nj, Nq); Kokkos::parallel_for( thread_policy, [&](const int i, const int j, const int q) { - const int q_point = - i + j * n_q_points_1d + q * n_q_points_1d * n_q_points_1d; - - // This loop simply multiplies the shape function at the quadrature - // point by the value finite element coefficient. - // FIXME check why using parallel_reduce ThreadVector is slower - const int base_shape = dof_to_quad ? q : q * n_q_points_1d; - const int stride_shape = dof_to_quad ? n_q_points_1d : 1; + const int base_shape = contract_over_rows ? q : q * n_columns; + const int stride_shape = contract_over_rows ? n_columns : 1; + const int base_in = - (direction == 0 ? - (n_q_points_1d * (i + n_q_points_1d * j)) : - (direction == 1 ? (i + n_q_points_1d * n_q_points_1d * j) : - (i + n_q_points_1d * j))); - const int stride_in = Utilities::pow(n_q_points_1d, direction); - Number sum = shape_data[base_shape] * in(base_in); - for (int k = 1; k < n_q_points_1d; ++k) - { - sum += shape_data[base_shape + k * stride_shape] * - in(base_in + k * stride_in); - } - t[q_point] = sum; - }); + (direction == 0 ? (i * Nj + j) * Nk : + (direction == 1 ? i + j * Ni * Nk : i * Nj + j)); + const int stride_in = Utilities::pow(n_columns, direction); - if (in_place) - team_member.team_barrier(); + Number sum = shape_data(base_shape) * in(base_in); + for (int k = 1; k < Nk; ++k) + sum += shape_data(base_shape + k * stride_shape) * + in(base_in + k * stride_in); - Kokkos::parallel_for( - thread_policy, [&](const int i, const int j, const int q) { - const int q_point = - i + j * n_q_points_1d + q * n_q_points_1d * n_q_points_1d; - const unsigned int destination_idx = - (direction == 0) ? (q + n_q_points_1d * (i + n_q_points_1d * j)) : - (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) : - (i + n_q_points_1d * (j + n_q_points_1d * q)); - - if (add) - Kokkos::atomic_add(&out(destination_idx), t[q_point]); + const int index_out = + (direction == 0 ? (i * Nj + j) * Nq + q : + (direction == 1 ? i + (j * Nq + q) * Ni : + (i + q * Ni) * Nj + j)); + + if constexpr (add) + Kokkos::atomic_add(&out(index_out), sum); else - out(destination_idx) = t[q_point]; + out(index_out) = sum; }); + + team_member.team_barrier(); } #endif - /** - * Helper function for values() and gradients(). - */ template DEAL_II_HOST_DEVICE void @@ -259,84 +306,68 @@ namespace Portable { #if KOKKOS_VERSION >= 40000 if constexpr (dim == 1) - apply_1d( + apply_1d( team_member, shape_data, in, out); if constexpr (dim == 2) - apply_2d( + apply_2d( team_member, shape_data, in, out); if constexpr (dim == 3) - apply_3d( + apply_3d( team_member, shape_data, in, out); #else - constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim); + // I: [0, m^{dim - direction - 1}) + // J: [0, n^direction) + constexpr int NI = Utilities::pow(n_rows, dim - direction - 1); + constexpr int NJ = Utilities::pow(n_columns, direction); + + constexpr int Nk = contract_over_rows ? n_rows : n_columns; + constexpr int Nq = contract_over_rows ? n_columns : n_rows; + + Assert(shape_data.size() == n_rows * n_columns, ExcInternalError()); + Assert(in.size() >= NI * NJ * Nk, ExcInternalError()); + Assert(out.size() >= NI * NJ * Nq, ExcInternalError()); + + constexpr int N = NI * NJ * Nq; + constexpr int stride = Utilities::pow(n_columns, direction); - Number t[n_q_points]; Kokkos::parallel_for( - Kokkos::TeamThreadRange(team_member, n_q_points), - [&](const int &q_point) { - const unsigned int i = (dim == 1) ? 0 : q_point % n_q_points_1d; - const unsigned int j = - (dim == 3) ? (q_point / n_q_points_1d) % n_q_points_1d : 0; - const unsigned int q = - (dim == 1) ? q_point : - (dim == 2) ? (q_point / n_q_points_1d) % n_q_points_1d : - q_point / (n_q_points_1d * n_q_points_1d); - - // This loop simply multiplies the shape function at the quadrature - // point by the value finite element coefficient. - const int stride_shape = dof_to_quad ? n_q_points_1d : 1; - const int stride = Utilities::pow(n_q_points_1d, direction); - const int base_shape = dof_to_quad ? q : (q * n_q_points_1d); - const int base = - (direction == 0) ? (n_q_points_1d * (i + n_q_points_1d * j)) : - (direction == 1) ? (i + n_q_points_1d * (n_q_points_1d * j)) : - (i + n_q_points_1d * j); - Number sum = - shape_data[base_shape] * (in_place ? out(base) : in(base)); - for (int k = 1; k < n_q_points_1d; ++k) + Kokkos::TeamThreadRange(team_member, N), [&](const int index_out) { + // index_in = (I Nk + k) n^direction + J + // index_out = (I Nq + q) n^direction + J + const int q = (index_out / stride) % Nq; + const int I = (index_out / stride) / Nq; + const int J = index_out % stride; + + const int base_shape = contract_over_rows ? q : q * n_columns; + const int stride_shape = contract_over_rows ? n_columns : 1; + const int base_in = I * Nk * stride + J; + + Number sum = shape_data(base_shape) * in(base_in); + for (int k = 1; k < Nk; ++k) { - sum += - shape_data[base_shape + k * stride_shape] * - (in_place ? out(base + k * stride) : in(base + k * stride)); + const int index_in = (I * Nk + k) * stride + J; + sum += shape_data(base_shape + k * stride_shape) * in(index_in); } - t[q_point] = sum; - }); - - if (in_place) - team_member.team_barrier(); - Kokkos::parallel_for( - Kokkos::TeamThreadRange(team_member, n_q_points), - [&](const int &q_point) { - const unsigned int i = (dim == 1) ? 0 : q_point % n_q_points_1d; - const unsigned int j = - (dim == 3) ? (q_point / n_q_points_1d) % n_q_points_1d : 0; - const unsigned int q = - (dim == 1) ? q_point : - (dim == 2) ? (q_point / n_q_points_1d) % n_q_points_1d : - q_point / (n_q_points_1d * n_q_points_1d); - - const unsigned int destination_idx = - (direction == 0) ? (q + n_q_points_1d * (i + n_q_points_1d * j)) : - (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) : - (i + n_q_points_1d * (j + n_q_points_1d * q)); - - if (add) - Kokkos::atomic_add(&out(destination_idx), t[q_point]); + if constexpr (add) + Kokkos::atomic_add(&out(index_out), sum); else - out(destination_idx) = t[q_point]; + out(index_out) = sum; }); + + team_member.team_barrier(); #endif } + /** * Generic evaluator framework. */ template struct EvaluatorTensorProduct {}; @@ -347,17 +378,22 @@ namespace Portable * Internal evaluator for 1d-3d shape function using the tensor product form * of the basis functions. */ - template + template struct EvaluatorTensorProduct { public: using TeamHandle = Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space>::member_type; + using SharedView = Kokkos::View>; + DEAL_II_HOST_DEVICE EvaluatorTensorProduct( const TeamHandle &team_member, @@ -365,53 +401,8 @@ namespace Portable Kokkos::View shape_gradients, Kokkos::View - co_shape_gradients); - - /** - * Evaluate the finite element function at the quadrature points. - */ - template - DEAL_II_HOST_DEVICE void - evaluate_values(ViewType u); - - /** - * Evaluate the gradients of the finite element function at the quadrature - * points. - */ - template - DEAL_II_HOST_DEVICE void - evaluate_gradients(const ViewTypeIn u, ViewTypeOut grad_u); - - /** - * Evaluate the values and the gradients of the finite element function at - * the quadrature points. - */ - template - DEAL_II_HOST_DEVICE void - evaluate_values_and_gradients(ViewType1 u, ViewType2 grad_u); - - /** - * Helper function for integrate(). Integrate the finite element function. - */ - template - DEAL_II_HOST_DEVICE void - integrate_values(ViewType u); - - /** - * Helper function for integrate(). Integrate the gradients of the finite - * element function. - */ - template - DEAL_II_HOST_DEVICE void - integrate_gradients(ViewType1 u, ViewType2 grad_u); - - /** - * Helper function for integrate(). Integrate the values and the gradients - * of the finite element function. - */ - template - DEAL_II_HOST_DEVICE void - integrate_values_and_gradients(ViewType1 u, ViewType2 grad_u); + co_shape_gradients, + SharedView temp); /** * Evaluate/integrate the values of a finite element function at the @@ -439,7 +430,6 @@ namespace Portable DEAL_II_HOST_DEVICE void gradients(const ViewTypeIn in, ViewTypeOut out) const; - public: /** * Evaluate the gradient of a finite element function at the quadrature * points for a given @p direction for collocation methods. @@ -474,33 +464,36 @@ namespace Portable */ Kokkos::View co_shape_gradients; + + /** + * Temporary storage for in-place evaluations. + */ + SharedView temp; }; - template + template DEAL_II_HOST_DEVICE - EvaluatorTensorProduct:: + EvaluatorTensorProduct:: EvaluatorTensorProduct( const TeamHandle &team_member, Kokkos::View shape_values, Kokkos::View shape_gradients, Kokkos::View - co_shape_gradients) + co_shape_gradients, + SharedView temp) : team_member(team_member) , shape_values(shape_values) , shape_gradients(shape_gradients) , co_shape_gradients(co_shape_gradients) + , temp(temp) {} - template + template template DEAL_II_HOST_DEVICE void - EvaluatorTensorProduct::values(const ViewTypeIn in, - ViewTypeOut out) const + EvaluatorTensorProduct:: + values(const ViewTypeIn in, ViewTypeOut out) const { - apply( - team_member, shape_values, in, out); + if constexpr (in_place) + { + apply( + team_member, shape_values, in, temp); + + populate_view(team_member, out, temp, out.extent(0)); + } + else + apply( + team_member, shape_values, in, out); } - template + template template DEAL_II_HOST_DEVICE void - EvaluatorTensorProduct::gradients(const ViewTypeIn in, - ViewTypeOut out) const + EvaluatorTensorProduct:: + gradients(const ViewTypeIn in, ViewTypeOut out) const { - apply( - team_member, shape_gradients, in, out); + if constexpr (in_place) + { + apply( + team_member, shape_gradients, in, temp); + + populate_view(team_member, out, temp, out.extent(0)); + } + else + apply( + team_member, shape_gradients, in, out); } - template + template template DEAL_II_HOST_DEVICE void - EvaluatorTensorProduct::co_gradients(const ViewTypeIn in, - ViewTypeOut out) const - { - apply( - team_member, co_shape_gradients, in, out); - } - - - - template - template - DEAL_II_HOST_DEVICE inline void - EvaluatorTensorProduct::evaluate_values(ViewType u) - { - if constexpr (dim == 1) - values<0, true, false, true>(u, u); - else if constexpr (dim == 2) - { - values<0, true, false, true>(u, u); - team_member.team_barrier(); - values<1, true, false, true>(u, u); - } - else if constexpr (dim == 3) - { - values<0, true, false, true>(u, u); - team_member.team_barrier(); - values<1, true, false, true>(u, u); - team_member.team_barrier(); - values<2, true, false, true>(u, u); - } - else - Kokkos::abort("dim must not exceed 3!"); - } - - - - template - template - DEAL_II_HOST_DEVICE inline void - EvaluatorTensorProduct::integrate_values(ViewType u) + EvaluatorTensorProduct:: + co_gradients(const ViewTypeIn in, ViewTypeOut out) const { - if constexpr (dim == 1) - values<0, false, false, true>(u, u); - else if constexpr (dim == 2) - { - values<0, false, false, true>(u, u); - team_member.team_barrier(); - values<1, false, false, true>(u, u); - } - else if constexpr (dim == 3) - { - values<0, false, false, true>(u, u); - team_member.team_barrier(); - values<1, false, false, true>(u, u); - team_member.team_barrier(); - values<2, false, false, true>(u, u); - } - else - Kokkos::abort("dim must not exceed 3!"); - } - - - - template - template - DEAL_II_HOST_DEVICE inline void - EvaluatorTensorProduct::evaluate_gradients(const ViewTypeIn u, - ViewTypeOut grad_u) - { - if constexpr (dim == 1) - { - gradients<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); - } - else if constexpr (dim == 2) - { - gradients<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); - values<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 1)); - - team_member.team_barrier(); - - values<1, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0), - Kokkos::subview(grad_u, Kokkos::ALL, 0)); - gradients<1, true, false, true>( - Kokkos::subview(grad_u, Kokkos::ALL, 1), - Kokkos::subview(grad_u, Kokkos::ALL, 1)); - } - else if constexpr (dim == 3) - { - gradients<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); - values<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 1)); - values<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 2)); - - team_member.team_barrier(); - - values<1, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0), - Kokkos::subview(grad_u, Kokkos::ALL, 0)); - gradients<1, true, false, true>( - Kokkos::subview(grad_u, Kokkos::ALL, 1), - Kokkos::subview(grad_u, Kokkos::ALL, 1)); - values<1, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 2), - Kokkos::subview(grad_u, Kokkos::ALL, 2)); - - team_member.team_barrier(); - - values<2, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0), - Kokkos::subview(grad_u, Kokkos::ALL, 0)); - values<2, true, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 1), - Kokkos::subview(grad_u, Kokkos::ALL, 1)); - gradients<2, true, false, true>( - Kokkos::subview(grad_u, Kokkos::ALL, 2), - Kokkos::subview(grad_u, Kokkos::ALL, 2)); - } - else - Kokkos::abort("dim must not exceed 3!"); - } - - - - template - template - DEAL_II_HOST_DEVICE inline void - EvaluatorTensorProduct::evaluate_values_and_gradients(ViewType1 u, - ViewType2 - grad_u) - { - if constexpr (dim == 1) - { - values<0, true, false, true>(u, u); - team_member.team_barrier(); - - co_gradients<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); - } - else if constexpr (dim == 2) - { - values<0, true, false, true>(u, u); - team_member.team_barrier(); - values<1, true, false, true>(u, u); - team_member.team_barrier(); - - co_gradients<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); - co_gradients<1, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 1)); - } - else if constexpr (dim == 3) - { - values<0, true, false, true>(u, u); - team_member.team_barrier(); - values<1, true, false, true>(u, u); - team_member.team_barrier(); - values<2, true, false, true>(u, u); - team_member.team_barrier(); - - co_gradients<0, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 0)); - co_gradients<1, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 1)); - co_gradients<2, true, false, false>( - u, Kokkos::subview(grad_u, Kokkos::ALL, 2)); - } - else - Kokkos::abort("dim must not exceed 3!"); - } - - - - template - template - DEAL_II_HOST_DEVICE inline void - EvaluatorTensorProduct::integrate_gradients(ViewType1 u, - ViewType2 grad_u) - { - if constexpr (dim == 1) - { - gradients<0, false, add, false>( - Kokkos::subview(grad_u, Kokkos::ALL, dim), u); - } - else if constexpr (dim == 2) - { - gradients<0, false, false, true>( - Kokkos::subview(grad_u, Kokkos::ALL, 0), - Kokkos::subview(grad_u, Kokkos::ALL, 0)); - values<0, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 1), - Kokkos::subview(grad_u, - Kokkos::ALL, - 1)); - - team_member.team_barrier(); - - values<1, false, add, false>(Kokkos::subview(grad_u, Kokkos::ALL, 0), - u); - team_member.team_barrier(); - gradients<1, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 1), u); - } - else if constexpr (dim == 3) - { - gradients<0, false, false, true>( - Kokkos::subview(grad_u, Kokkos::ALL, 0), - Kokkos::subview(grad_u, Kokkos::ALL, 0)); - values<0, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 1), - Kokkos::subview(grad_u, - Kokkos::ALL, - 1)); - values<0, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 2), - Kokkos::subview(grad_u, - Kokkos::ALL, - 2)); - - team_member.team_barrier(); - - values<1, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 0), - Kokkos::subview(grad_u, - Kokkos::ALL, - 0)); - gradients<1, false, false, true>( - Kokkos::subview(grad_u, Kokkos::ALL, 1), - Kokkos::subview(grad_u, Kokkos::ALL, 1)); - values<1, false, false, true>(Kokkos::subview(grad_u, Kokkos::ALL, 2), - Kokkos::subview(grad_u, - Kokkos::ALL, - 2)); - - team_member.team_barrier(); - - values<2, false, add, false>(Kokkos::subview(grad_u, Kokkos::ALL, 0), - u); - team_member.team_barrier(); - values<2, false, true, false>(Kokkos::subview(grad_u, Kokkos::ALL, 1), - u); - team_member.team_barrier(); - gradients<2, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 2), u); - } - else - Kokkos::abort("dim must not exceed 3!"); - } - - - - template - template - DEAL_II_HOST_DEVICE inline void - EvaluatorTensorProduct::integrate_values_and_gradients(ViewType1 u, - ViewType2 - grad_u) - { - if constexpr (dim == 1) - { - co_gradients<0, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 0), u); - team_member.team_barrier(); - - values<0, false, false, true>(u, u); - } - else if constexpr (dim == 2) - { - co_gradients<1, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 1), u); - team_member.team_barrier(); - co_gradients<0, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 0), u); - team_member.team_barrier(); - - values<1, false, false, true>(u, u); - team_member.team_barrier(); - values<0, false, false, true>(u, u); - team_member.team_barrier(); - } - else if constexpr (dim == 3) + if constexpr (in_place) { - co_gradients<2, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 2), u); - team_member.team_barrier(); - co_gradients<1, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 1), u); - team_member.team_barrier(); - co_gradients<0, false, true, false>( - Kokkos::subview(grad_u, Kokkos::ALL, 0), u); - team_member.team_barrier(); - - values<2, false, false, true>(u, u); - team_member.team_barrier(); - values<1, false, false, true>(u, u); - team_member.team_barrier(); - values<0, false, false, true>(u, u); - team_member.team_barrier(); + apply(team_member, co_shape_gradients, in, temp); + + populate_view(team_member, out, temp, out.extent(0)); } else - Kokkos::abort("dim must not exceed 3!"); + apply( + team_member, co_shape_gradients, in, out); } } // namespace internal } // namespace Portable diff --git a/tests/matrix_free_kokkos/evaluate_1d_shape.cc b/tests/matrix_free_kokkos/evaluate_1d_shape.cc index 7daa538697..037bd42109 100644 --- a/tests/matrix_free_kokkos/evaluate_1d_shape.cc +++ b/tests/matrix_free_kokkos/evaluate_1d_shape.cc @@ -44,13 +44,23 @@ evaluate_tensor_product( Kokkos::View dst, Kokkos::View src) { + Kokkos::View< + double *, + MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, + Kokkos::MemoryTraits> + dummy_scratch(team_member.team_shmem(), 0); + Portable::internal::EvaluatorTensorProduct< Portable::internal::evaluate_general, 1, - M - 1, + M, N, double> - evaluator(team_member, shape_values, shape_gradients, co_shape_gradients); + evaluator(team_member, + shape_values, + shape_gradients, + co_shape_gradients, + dummy_scratch); if (type == 0) evaluator.template values<0, dof_to_quad, add, false>(src, dst); @@ -197,19 +207,19 @@ main() deallog.push("values"); test<4, 4, 0, false>(); test<3, 3, 0, false>(); + test<3, 4, 0, false>(); + test<3, 5, 0, false>(); // Not supported right now // test<4,3,0,false>(); - // test<3,4,0,false>(); - // test<3,5,0,false>(); deallog.pop(); deallog.push("gradients"); test<4, 4, 1, false>(); test<3, 3, 1, false>(); + test<3, 4, 1, false>(); + test<3, 5, 1, false>(); // Not supported right now // test<4,3,1,false>(); - // test<3,4,1,false>(); - // test<3,5,1,false>(); deallog.pop(); deallog.push("add"); @@ -217,19 +227,19 @@ main() deallog.push("values"); test<4, 4, 0, true>(); test<3, 3, 0, true>(); + test<3, 4, 0, true>(); + test<3, 5, 0, true>(); // Not supported right now // test<4,3,0,true>(); - // test<3,4,0,true>(); - // test<3,5,0,true>(); deallog.pop(); deallog.push("gradients"); test<4, 4, 1, true>(); test<3, 3, 1, true>(); + test<3, 4, 1, true>(); + test<3, 5, 1, true>(); // Not supported right now // test<4,3,1,true>(); - // test<3,4,1,true>(); - // test<3,5,1,true>(); deallog.pop(); deallog.pop(); diff --git a/tests/matrix_free_kokkos/evaluate_1d_shape.output b/tests/matrix_free_kokkos/evaluate_1d_shape.output index 5aa83953a8..b18c6ceb08 100644 --- a/tests/matrix_free_kokkos/evaluate_1d_shape.output +++ b/tests/matrix_free_kokkos/evaluate_1d_shape.output @@ -5,21 +5,45 @@ DEAL:values::Errors transpose: 0 0 0 0 DEAL:values::Test 3 x 3 DEAL:values::Errors no transpose: 0 0 0 DEAL:values::Errors transpose: 0 0 0 +DEAL:values::Test 3 x 4 +DEAL:values::Errors no transpose: 0 0 0 +DEAL:values::Errors transpose: 0 0 0 0 +DEAL:values::Test 3 x 5 +DEAL:values::Errors no transpose: 0 0 0 +DEAL:values::Errors transpose: 0 0 0 0 0 DEAL:gradients::Test 4 x 4 DEAL:gradients::Errors no transpose: 0 0 0 0 DEAL:gradients::Errors transpose: 0 0 0 0 DEAL:gradients::Test 3 x 3 DEAL:gradients::Errors no transpose: 0 0 0 DEAL:gradients::Errors transpose: 0 0 0 +DEAL:gradients::Test 3 x 4 +DEAL:gradients::Errors no transpose: 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 +DEAL:gradients::Test 3 x 5 +DEAL:gradients::Errors no transpose: 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 0 DEAL:add:values::Test 4 x 4 DEAL:add:values::Errors no transpose: 0 0 0 0 DEAL:add:values::Errors transpose: 0 0 0 0 DEAL:add:values::Test 3 x 3 DEAL:add:values::Errors no transpose: 0 0 0 DEAL:add:values::Errors transpose: 0 0 0 +DEAL:add:values::Test 3 x 4 +DEAL:add:values::Errors no transpose: 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 +DEAL:add:values::Test 3 x 5 +DEAL:add:values::Errors no transpose: 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 0 DEAL:add:gradients::Test 4 x 4 DEAL:add:gradients::Errors no transpose: 0 0 0 0 DEAL:add:gradients::Errors transpose: 0 0 0 0 DEAL:add:gradients::Test 3 x 3 DEAL:add:gradients::Errors no transpose: 0 0 0 DEAL:add:gradients::Errors transpose: 0 0 0 +DEAL:add:gradients::Test 3 x 4 +DEAL:add:gradients::Errors no transpose: 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 +DEAL:add:gradients::Test 3 x 5 +DEAL:add:gradients::Errors no transpose: 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 0 diff --git a/tests/matrix_free_kokkos/evaluate_2d_shape.cc b/tests/matrix_free_kokkos/evaluate_2d_shape.cc index a2f51058e4..e82a4fea4f 100644 --- a/tests/matrix_free_kokkos/evaluate_2d_shape.cc +++ b/tests/matrix_free_kokkos/evaluate_2d_shape.cc @@ -41,27 +41,41 @@ evaluate_tensor_product( Kokkos::View shape_gradients, Kokkos::View co_shape_gradients, Kokkos::View dst, - Kokkos::View src) + Kokkos::View src, + Kokkos::View tmp) { + Kokkos::View< + double *, + MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, + Kokkos::MemoryTraits> + dummy_scratch(team_member.team_shmem(), 0); + Portable::internal::EvaluatorTensorProduct< Portable::internal::evaluate_general, 2, - M - 1, + M, N, double> - evaluator(team_member, shape_values, shape_gradients, co_shape_gradients); + evaluator(team_member, + shape_values, + shape_gradients, + co_shape_gradients, + dummy_scratch); + + constexpr int d0 = dof_to_quad ? 0 : 1; + constexpr int d1 = dof_to_quad ? 1 : 0; if (type == 0) { - evaluator.template values<0, dof_to_quad, false, false>(src, src); + evaluator.template values(src, tmp); team_member.team_barrier(); - evaluator.template values<1, dof_to_quad, add, false>(src, dst); + evaluator.template values(tmp, dst); } if (type == 1) { - evaluator.template gradients<0, dof_to_quad, false, false>(src, src); + evaluator.template gradients(src, tmp); team_member.team_barrier(); - evaluator.template gradients<1, dof_to_quad, add, false>(src, dst); + evaluator.template gradients(tmp, dst); } } @@ -93,11 +107,14 @@ test() constexpr int M_2d = M * M; constexpr int N_2d = N * N; + constexpr int MN = M * N; LinearAlgebra::ReadWriteVector x_ref(N_2d), y_ref(M_2d); Kokkos::View x_dev( Kokkos::view_alloc("x_dev", Kokkos::WithoutInitializing)); Kokkos::View y_dev( Kokkos::view_alloc("y_dev", Kokkos::WithoutInitializing)); + Kokkos::View tmp_dev( + Kokkos::view_alloc("tmp_dev", Kokkos::WithoutInitializing)); auto x_host = Kokkos::create_mirror_view(x_dev); auto y_host = Kokkos::create_mirror_view(y_dev); @@ -160,7 +177,8 @@ test() shape_gradients, co_shape_gradients, y_dev, - x_dev); + x_dev, + tmp_dev); }); // Check the results on the host @@ -197,7 +215,8 @@ test() shape_gradients, co_shape_gradients, x_dev, - y_dev); + y_dev, + tmp_dev); }); // Check the results on the host @@ -219,11 +238,15 @@ main() deallog.push("values"); test<4, 4, 0, false>(); test<3, 3, 0, false>(); + test<3, 4, 0, false>(); + test<3, 5, 0, false>(); deallog.pop(); deallog.push("gradients"); test<4, 4, 1, false>(); test<3, 3, 1, false>(); + test<3, 4, 1, false>(); + test<3, 5, 1, false>(); deallog.pop(); deallog.push("add"); @@ -231,11 +254,15 @@ main() deallog.push("values"); test<4, 4, 0, true>(); test<3, 3, 0, true>(); + test<3, 4, 0, true>(); + test<3, 5, 0, true>(); deallog.pop(); deallog.push("gradients"); test<4, 4, 1, true>(); test<3, 3, 1, true>(); + test<3, 4, 1, true>(); + test<3, 5, 1, true>(); deallog.pop(); deallog.pop(); diff --git a/tests/matrix_free_kokkos/evaluate_2d_shape.output b/tests/matrix_free_kokkos/evaluate_2d_shape.output index 312a78a60e..432869ab0c 100644 --- a/tests/matrix_free_kokkos/evaluate_2d_shape.output +++ b/tests/matrix_free_kokkos/evaluate_2d_shape.output @@ -5,21 +5,45 @@ DEAL:values::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:values::Test 3 x 3 DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 DEAL:values::Errors transpose: 0 0 0 0 0 0 0 0 0 +DEAL:values::Test 3 x 4 +DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:values::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:values::Test 3 x 5 +DEAL:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:values::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:gradients::Test 4 x 4 DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:gradients::Test 3 x 3 -DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 -DEAL:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 +DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 +DEAL:gradients::Test 3 x 4 +DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:gradients::Test 3 x 5 +DEAL:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:add:values::Test 4 x 4 -DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:add:values::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:add:values::Test 3 x 3 -DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 -DEAL:add:values::Errors transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Test 3 x 4 +DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Test 3 x 5 +DEAL:add:values::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:values::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:add:gradients::Test 4 x 4 -DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -DEAL:add:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 DEAL:add:gradients::Test 3 x 3 -DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 -DEAL:add:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Test 3 x 4 +DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Test 3 x 5 +DEAL:add:gradients::Errors no transpose: 0 0 0 0 0 0 0 0 0 +DEAL:add:gradients::Errors transpose: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/tests/tests.h b/tests/tests.h index 49ff8f7c7c..b9b5fa3678 100644 --- a/tests/tests.h +++ b/tests/tests.h @@ -72,9 +72,13 @@ struct DisableWindowsDebugRuntimeDialog // Redefine Assert as AssertThrow to make sure that the code is tested similarly // in Release mode and in Debug mode. clang-format makes sure that this file is // included after all regular header files but before all the other local header -// files. -#undef Assert -#define Assert AssertThrow +// files. The redefinition is not applied when Kokkos GPU backend is enabled, +// because AssertThrow is not defined for device codes. +#if !defined(KOKKOS_ENABLE_CUDA) && !defined(KOKKOS_ENABLE_HIP) && \ + !defined(KOKKOS_ENABLE_SYCL) +# undef Assert +# define Assert AssertThrow +#endif // implicitly use the deal.II namespace everywhere, without us having to say // so in each and every testcase -- 2.39.5