From c05614646d785160d6ae3f1916ec6133a9f7e593 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Tue, 13 Aug 2024 16:46:05 +0200 Subject: [PATCH] Support multiple components in Portable::MatrixFree. --- doc/news/changes/minor/20240814Munch | 4 + .../matrix_free/portable_fe_evaluation.h | 259 ++++++++---- .../portable_hanging_nodes_internal.h | 33 +- .../matrix_free/portable_matrix_free.h | 33 +- .../portable_matrix_free.templates.h | 39 +- ...trix_vector_host_device_multi_component.cc | 383 ++++++++++++++++++ ..._vector_host_device_multi_component.output | 20 + 7 files changed, 651 insertions(+), 120 deletions(-) create mode 100644 doc/news/changes/minor/20240814Munch create mode 100644 tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.cc create mode 100644 tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.output diff --git a/doc/news/changes/minor/20240814Munch b/doc/news/changes/minor/20240814Munch new file mode 100644 index 0000000000..59da529b15 --- /dev/null +++ b/doc/news/changes/minor/20240814Munch @@ -0,0 +1,4 @@ +Improved: Portable::MatrixFree now supports FESystem(FE_Q, n_components) +with arbitrary number of components. +
+(Peter Munch, Daniel Arndt, 2024/08/14) diff --git a/include/deal.II/matrix_free/portable_fe_evaluation.h b/include/deal.II/matrix_free/portable_fe_evaluation.h index f94adf548d..d2db2377ca 100644 --- a/include/deal.II/matrix_free/portable_fe_evaluation.h +++ b/include/deal.II/matrix_free/portable_fe_evaluation.h @@ -72,12 +72,19 @@ namespace Portable /** * An alias for scalar quantities. */ - using value_type = Number; + using value_type = std::conditional_t<(n_components_ == 1), + Number, + Tensor<1, n_components_, Number>>; /** * An alias for vectorial quantities. */ - using gradient_type = Tensor<1, dim, Number>; + using gradient_type = std::conditional_t< + n_components_ == 1, + Tensor<1, dim, Number>, + std::conditional_t, + Tensor<1, n_components_, Tensor<1, dim, Number>>>>; /** * An alias to kernel specific information. @@ -265,22 +272,24 @@ namespace Portable FEEvaluation:: read_dof_values(const Number *src) { - static_assert(n_components_ == 1, "This function only supports FE with one \ - components"); // Populate the scratch memory - Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member, - n_q_points), - [&](const int &i) { - shared_data->values(i) = - src[data->local_to_global(cell_id, i)]; - }); + 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 * c)]; + }); shared_data->team_member.team_barrier(); - internal::resolve_hanging_nodes( - shared_data->team_member, - data->constraint_weights, - data->constraint_mask(cell_id), - shared_data->values); + for (unsigned int c = 0; c < n_components_; ++c) + { + internal::resolve_hanging_nodes( + shared_data->team_member, + data->constraint_weights, + data->constraint_mask(cell_id), + Kokkos::subview(shared_data->values, Kokkos::ALL, c)); + } } @@ -294,31 +303,35 @@ namespace Portable FEEvaluation:: distribute_local_to_global(Number *dst) const { - static_assert(n_components_ == 1, "This function only supports FE with one \ - components"); - - internal::resolve_hanging_nodes( - shared_data->team_member, - data->constraint_weights, - data->constraint_mask(cell_id), - shared_data->values); + for (unsigned int c = 0; c < n_components_; ++c) + { + internal::resolve_hanging_nodes( + shared_data->team_member, + data->constraint_weights, + data->constraint_mask(cell_id), + Kokkos::subview(shared_data->values, Kokkos::ALL, c)); + } if (data->use_coloring) { - Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member, - n_q_points), - [&](const int &i) { - dst[data->local_to_global(cell_id, i)] += - shared_data->values(i); - }); + Kokkos::parallel_for( + Kokkos::TeamThreadRange(shared_data->team_member, n_q_points), + [&](const int &i) { + for (unsigned int c = 0; c < n_components_; ++c) + dst[data->local_to_global(cell_id, + i + tensor_dofs_per_cell * c)] += + shared_data->values(i, c); + }); } else { Kokkos::parallel_for( Kokkos::TeamThreadRange(shared_data->team_member, n_q_points), [&](const int &i) { - Kokkos::atomic_add(&dst[data->local_to_global(cell_id, i)], - shared_data->values(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 * c)], + shared_data->values(i, c)); }); } } @@ -347,23 +360,31 @@ namespace Portable data->shape_gradients, data->co_shape_gradients); - if ((evaluate_flag & EvaluationFlags::values) && - (evaluate_flag & EvaluationFlags::gradients)) - { - evaluator_tensor_product.evaluate_values_and_gradients( - shared_data->values, shared_data->gradients); - shared_data->team_member.team_barrier(); - } - else if (evaluate_flag & EvaluationFlags::gradients) - { - evaluator_tensor_product.evaluate_gradients(shared_data->values, - shared_data->gradients); - shared_data->team_member.team_barrier(); - } - else if (evaluate_flag & EvaluationFlags::values) + for (unsigned int c = 0; c < n_components_; ++c) { - evaluator_tensor_product.evaluate_values(shared_data->values); - shared_data->team_member.team_barrier(); + 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(); + } } } @@ -406,22 +427,31 @@ namespace Portable data->shape_gradients, data->co_shape_gradients); - if ((integration_flag & EvaluationFlags::values) && - (integration_flag & EvaluationFlags::gradients)) - { - evaluator_tensor_product.integrate_values_and_gradients( - shared_data->values, shared_data->gradients); - } - else if (integration_flag & EvaluationFlags::values) - { - evaluator_tensor_product.integrate_values(shared_data->values); - shared_data->team_member.team_barrier(); - } - else if (integration_flag & EvaluationFlags::gradients) + + for (unsigned int c = 0; c < n_components_; ++c) { - evaluator_tensor_product.template integrate_gradients( - shared_data->values, shared_data->gradients); - shared_data->team_member.team_barrier(); + 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(); + } } } @@ -457,7 +487,17 @@ namespace Portable FEEvaluation::get_value( int q_point) const { - return shared_data->values(q_point); + if constexpr (n_components_ == 1) + { + return shared_data->values(q_point, 0); + } + else + { + value_type result; + for (unsigned int c = 0; c < n_components; ++c) + result[c] = shared_data->values(q_point, c); + return result; + } } @@ -475,7 +515,17 @@ namespace Portable FEEvaluation:: get_dof_value(int q_point) const { - return shared_data->values(q_point); + if constexpr (n_components_ == 1) + { + return shared_data->values(q_point, 0); + } + else + { + value_type result; + for (unsigned int c = 0; c < n_components; ++c) + result[c] = shared_data->values(q_point, c); + return result; + } } @@ -489,7 +539,16 @@ namespace Portable FEEvaluation:: submit_value(const value_type &val_in, int q_point) { - shared_data->values(q_point) = val_in * data->JxW(cell_id, q_point); + if constexpr (n_components_ == 1) + { + shared_data->values(q_point, 0) = val_in * data->JxW(cell_id, q_point); + } + else + { + for (unsigned int c = 0; c < n_components; ++c) + shared_data->values(q_point, c) = + val_in[c] * data->JxW(cell_id, q_point); + } } @@ -503,7 +562,15 @@ namespace Portable FEEvaluation:: submit_dof_value(const value_type &val_in, int q_point) { - shared_data->values(q_point) = val_in; + if constexpr (n_components_ == 1) + { + shared_data->values(q_point, 0) = val_in; + } + else + { + for (unsigned int c = 0; c < n_components; ++c) + shared_data->values(q_point, c) = val_in[c]; + } } @@ -521,17 +588,30 @@ namespace Portable FEEvaluation:: get_gradient(int q_point) const { - static_assert(n_components_ == 1, "This function only supports FE with one \ - components"); - gradient_type grad; - for (unsigned int d_1 = 0; d_1 < dim; ++d_1) + + if constexpr (n_components_ == 1) { - Number tmp = 0.; - for (unsigned int d_2 = 0; d_2 < dim; ++d_2) - tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) * - shared_data->gradients(q_point, d_2); - grad[d_1] = tmp; + for (unsigned int d_1 = 0; d_1 < dim; ++d_1) + { + Number tmp = 0.; + for (unsigned int d_2 = 0; d_2 < dim; ++d_2) + tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) * + shared_data->gradients(q_point, d_2, 0); + grad[d_1] = tmp; + } + } + else + { + for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int d_1 = 0; d_1 < dim; ++d_1) + { + Number tmp = 0.; + for (unsigned int d_2 = 0; d_2 < dim; ++d_2) + tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) * + shared_data->gradients(q_point, d_2, c); + grad[c][d_1] = tmp; + } } return grad; @@ -548,13 +628,30 @@ namespace Portable FEEvaluation:: submit_gradient(const gradient_type &grad_in, int q_point) { - for (unsigned int d_1 = 0; d_1 < dim; ++d_1) + if constexpr (n_components_ == 1) + { + for (unsigned int d_1 = 0; d_1 < dim; ++d_1) + { + Number tmp = 0.; + for (unsigned int d_2 = 0; d_2 < dim; ++d_2) + tmp += + data->inv_jacobian(cell_id, q_point, d_1, d_2) * grad_in[d_2]; + shared_data->gradients(q_point, d_1, 0) = + tmp * data->JxW(cell_id, q_point); + } + } + else { - Number tmp = 0.; - for (unsigned int d_2 = 0; d_2 < dim; ++d_2) - tmp += data->inv_jacobian(cell_id, q_point, d_1, d_2) * grad_in[d_2]; - shared_data->gradients(q_point, d_1) = - tmp * data->JxW(cell_id, q_point); + for (unsigned int c = 0; c < n_components; ++c) + for (unsigned int d_1 = 0; d_1 < dim; ++d_1) + { + Number tmp = 0.; + for (unsigned int d_2 = 0; d_2 < dim; ++d_2) + tmp += data->inv_jacobian(cell_id, q_point, d_1, d_2) * + grad_in[c][d_2]; + shared_data->gradients(q_point, d_1, c) = + tmp * data->JxW(cell_id, q_point); + } } } diff --git a/include/deal.II/matrix_free/portable_hanging_nodes_internal.h b/include/deal.II/matrix_free/portable_hanging_nodes_internal.h index 3139473b3f..d54dc72862 100644 --- a/include/deal.II/matrix_free/portable_hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/portable_hanging_nodes_internal.h @@ -121,7 +121,8 @@ namespace Portable template + typename Number, + typename ViewType> DEAL_II_HOST_DEVICE inline void interpolate_boundary_2d( const Kokkos::TeamPolicy< @@ -130,11 +131,8 @@ namespace Portable Kokkos::View constraint_weights, const dealii::internal::MatrixFreeFunctions::ConstraintKinds - &constraint_mask, - Kokkos::View> values) + &constraint_mask, + ViewType values) { constexpr unsigned int n_q_points_1d = fe_degree + 1; constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 2); @@ -242,7 +240,8 @@ namespace Portable template + typename Number, + typename ViewType> DEAL_II_HOST_DEVICE inline void interpolate_boundary_3d( const Kokkos::TeamPolicy< @@ -251,11 +250,8 @@ namespace Portable Kokkos::View constraint_weights, const dealii::internal::MatrixFreeFunctions::ConstraintKinds - constraint_mask, - Kokkos::View> values) + constraint_mask, + ViewType values) { constexpr unsigned int n_q_points_1d = fe_degree + 1; constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 3); @@ -410,7 +406,11 @@ namespace Portable * @cite ljungkvist2017matrix and in Section 3.4 of * @cite kronbichler2019multigrid. */ - template + template DEAL_II_HOST_DEVICE void resolve_hanging_nodes( const Kokkos::TeamPolicy< @@ -419,11 +419,8 @@ namespace Portable Kokkos::View constraint_weights, const dealii::internal::MatrixFreeFunctions::ConstraintKinds - constraint_mask, - Kokkos::View> values) + constraint_mask, + ViewType values) { if constexpr (dim == 2) { diff --git a/include/deal.II/matrix_free/portable_matrix_free.h b/include/deal.II/matrix_free/portable_matrix_free.h index 2b25b1ef3f..73551b2d31 100644 --- a/include/deal.II/matrix_free/portable_matrix_free.h +++ b/include/deal.II/matrix_free/portable_matrix_free.h @@ -217,6 +217,11 @@ namespace Portable */ unsigned int n_cells; + /** + * Number of components. + */ + unsigned int n_components; + /** * Length of the padding. */ @@ -504,6 +509,16 @@ namespace Portable */ unsigned int fe_degree; + /** + * Number of components. + */ + unsigned int n_components; + + /** + * Number of scalar degrees of freedom per cell. + */ + unsigned int scalar_dofs_per_cell; + /** * Number of degrees of freedom per cell. */ @@ -634,19 +649,19 @@ namespace Portable using TeamHandle = Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space>::member_type; - using SharedView1D = Kokkos::View< - Number *, + using SharedViewValues = Kokkos::View< + Number **, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits>; - using SharedView2D = Kokkos::View< - Number *[dim], + using SharedViewGradients = Kokkos::View< + Number ***, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits>; DEAL_II_HOST_DEVICE - SharedData(const TeamHandle &team_member, - const SharedView1D &values, - const SharedView2D &gradients) + SharedData(const TeamHandle &team_member, + const SharedViewValues &values, + const SharedViewGradients &gradients) : team_member(team_member) , values(values) , gradients(gradients) @@ -660,12 +675,12 @@ namespace Portable /** * Memory for dof and quad values. */ - SharedView1D values; + SharedViewValues values; /** * Memory for computed gradients in reference coordinate system. */ - SharedView2D gradients; + SharedViewGradients gradients; }; 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 4c8a3db6d1..ebbaeef14d 100644 --- a/include/deal.II/matrix_free/portable_matrix_free.templates.h +++ b/include/deal.II/matrix_free/portable_matrix_free.templates.h @@ -82,6 +82,8 @@ namespace Portable const std::vector &lexicographic_inv; std::vector lexicographic_dof_indices; const unsigned int fe_degree; + const unsigned int n_components; + const unsigned int scalar_dofs_per_cell; const unsigned int dofs_per_cell; const unsigned int q_points_per_cell; const UpdateFlags &update_flags; @@ -109,6 +111,8 @@ namespace Portable update_values | update_gradients | update_JxW_values) , lexicographic_inv(shape_info.lexicographic_numbering) , fe_degree(data->fe_degree) + , n_components(data->n_components) + , scalar_dofs_per_cell(data->scalar_dofs_per_cell) , dofs_per_cell(data->dofs_per_cell) , q_points_per_cell(data->q_points_per_cell) , update_flags(update_flags) @@ -180,7 +184,7 @@ namespace Portable Kokkos::view_alloc("JxW_" + std::to_string(color), Kokkos::WithoutInitializing), n_cells, - dofs_per_cell); + scalar_dofs_per_cell); if (update_flags & update_gradients) data->inv_jacobian[color] = @@ -188,7 +192,7 @@ namespace Portable Kokkos::view_alloc("inv_jacobian_" + std::to_string(color), Kokkos::WithoutInitializing), n_cells, - dofs_per_cell); + scalar_dofs_per_cell); // Initialize to zero, i.e., unconstrained cell data->constraint_mask[color] = @@ -354,13 +358,13 @@ namespace Portable { using TeamHandle = Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space>::member_type; - using SharedView1D = - Kokkos::View>; - using SharedView2D = - Kokkos::View>; @@ -386,8 +390,11 @@ namespace Portable size_t team_shmem_size(int /*team_size*/) const { - return SharedView1D::shmem_size(Functor::n_local_dofs) + - SharedView2D::shmem_size(Functor::n_local_dofs); + return SharedViewValues::shmem_size(Functor::n_local_dofs, + gpu_data.n_components) + + SharedViewGradients::shmem_size(Functor::n_local_dofs, + dim, + gpu_data.n_components); } @@ -396,8 +403,13 @@ namespace Portable operator()(const TeamHandle &team_member) const { // Get the scratch memory - SharedView1D values(team_member.team_shmem(), Functor::n_local_dofs); - SharedView2D gradients(team_member.team_shmem(), Functor::n_local_dofs); + SharedViewValues values(team_member.team_shmem(), + Functor::n_local_dofs, + gpu_data.n_components); + SharedViewGradients gradients(team_member.team_shmem(), + Functor::n_local_dofs, + dim, + gpu_data.n_components); SharedData shared_data(team_member, values, gradients); func(team_member.league_rank(), &gpu_data, &shared_data, src, dst); @@ -504,6 +516,7 @@ namespace Portable data_copy.co_shape_gradients = co_shape_gradients; data_copy.constraint_weights = constraint_weights; data_copy.n_cells = n_cells[color]; + data_copy.n_components = n_components; data_copy.padding_length = padding_length; data_copy.row_start = row_start[color]; data_copy.use_coloring = use_coloring; @@ -716,8 +729,10 @@ namespace Portable padding_length = 1 << static_cast( std::ceil(dim * std::log2(fe_degree + 1.))); - dofs_per_cell = fe.n_dofs_per_cell(); - q_points_per_cell = Utilities::fixed_power(n_q_points_1d); + dofs_per_cell = fe.n_dofs_per_cell(); + n_components = fe.n_components(); + scalar_dofs_per_cell = dofs_per_cell / n_components; + q_points_per_cell = Utilities::fixed_power(n_q_points_1d); ::dealii::internal::MatrixFreeFunctions::ShapeInfo shape_info(quad, fe); diff --git a/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.cc b/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.cc new file mode 100644 index 0000000000..10639ffba5 --- /dev/null +++ b/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.cc @@ -0,0 +1,383 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2019 - 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. +// +// ------------------------------------------------------------------------ + + + +// Test vector valued FEM for Portable::MatrixFree and Porable::FEEvaluation. + +#include + +#include +#include +#include + +#include + +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include "../tests.h" + +#include "Kokkos_Core.hpp" + +template +class LaplaceOperator; + +template +class LaplaceOperator +{ +public: + using VectorType = + LinearAlgebra::distributed::Vector; + + LaplaceOperator(const unsigned int version) + : version(version) + {} + + void + reinit(const Mapping &mapping, + const DoFHandler &dof_handler, + const AffineConstraints &constraints, + const Quadrature<1> &quadrature) + { + typename MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = update_gradients; + + matrix_free.reinit( + mapping, dof_handler, constraints, quadrature, additional_data); + } + + void + initialize_dof_vector(VectorType &vec) const + { + matrix_free.initialize_dof_vector(vec); + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + matrix_free.cell_loop(&LaplaceOperator::local_apply, this, dst, src, true); + } + +private: + void + local_apply(const MatrixFree &data, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const + { + FEEvaluation phi(data); + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values_plain(src); + + if (version == 0) + phi.evaluate(EvaluationFlags::values); + else if (version == 1) + phi.evaluate(EvaluationFlags::gradients); + else if (version == 2) + phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + + for (unsigned int q = 0; q < phi.n_q_points; ++q) + { + if (version == 0 || version == 2) + phi.submit_value(phi.get_value(q), q); + if (version == 1 || version == 2) + phi.submit_gradient(phi.get_gradient(q), q); + } + + if (version == 0) + phi.integrate(EvaluationFlags::values); + else if (version == 1) + phi.integrate(EvaluationFlags::gradients); + else if (version == 2) + phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + + phi.distribute_local_to_global(dst); + } + } + + const unsigned int version; + MatrixFree matrix_free; +}; + + + +template +class LaplaceOperatorQuad +{ +public: + DEAL_II_HOST_DEVICE + LaplaceOperatorQuad(const unsigned int version) + : version(version) + {} + + DEAL_II_HOST_DEVICE void + operator()( + Portable::FEEvaluation + *phi, + const int q) const + { + if (version == 0 || version == 2) + phi->submit_value(phi->get_value(q), q); + if (version == 1 || version == 2) + phi->submit_gradient(phi->get_gradient(q), q); + } + +private: + const unsigned int version; +}; + +template +class LaplaceOperatorLocal +{ +public: + LaplaceOperatorLocal(const unsigned int version) + : version(version) + {} + + DEAL_II_HOST_DEVICE void + operator()(const unsigned int cell, + const typename Portable::MatrixFree::Data *gpu_data, + Portable::SharedData *shared_data, + const Number *src, + Number *dst) const + { + (void)cell; // TODO? + + Portable::FEEvaluation + phi( + /*cell,*/ gpu_data, shared_data); + phi.read_dof_values(src); + + if (version == 0) + phi.evaluate(EvaluationFlags::values); + else if (version == 1) + phi.evaluate(EvaluationFlags::gradients); + else if (version == 2) + phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + + phi.apply_for_each_quad_point( + LaplaceOperatorQuad(version)); + + if (version == 0) + phi.integrate(EvaluationFlags::values); + else if (version == 1) + phi.integrate(EvaluationFlags::gradients); + else if (version == 2) + phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + + phi.distribute_local_to_global(dst); + } + static const unsigned int n_dofs_1d = fe_degree + 1; + static const unsigned int n_local_dofs = Utilities::pow(fe_degree + 1, dim); + static const unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim); + +private: + const unsigned int version; +}; + +template +class LaplaceOperator +{ +public: + using VectorType = + LinearAlgebra::distributed::Vector; + + LaplaceOperator(const unsigned int version) + : version(version) + {} + + void + reinit(const Mapping &mapping, + const DoFHandler &dof_handler, + const AffineConstraints &constraints, + const Quadrature<1> &quadrature) + { + typename Portable::MatrixFree::AdditionalData additional_data; + additional_data.mapping_update_flags = update_JxW_values | update_gradients; + + matrix_free.reinit( + mapping, dof_handler, constraints, quadrature, additional_data); + } + + void + initialize_dof_vector(VectorType &vec) const + { + matrix_free.initialize_dof_vector(vec); + } + + void + vmult(VectorType &dst, const VectorType &src) const + { + dst = 0.0; // TODO: annoying + LaplaceOperatorLocal local_operator( + version); + matrix_free.cell_loop(local_operator, src, dst); + matrix_free.copy_constrained_values(src, dst); // TODO: annoying + } + +private: + const unsigned int version; + Portable::MatrixFree matrix_free; +}; + + + +template +class AnalyticalFunction : public Function +{ +public: + AnalyticalFunction(const unsigned int n_components) + : Function(n_components) + {} + + virtual T + value(const Point &p, const unsigned int component = 0) const override + { + double temp = 0.0; + + for (unsigned int d = 0; d < dim; ++d) + temp += std::sin(p[d]); + + return temp * (1.0 + component); + } +}; + + + +template +void +run(const unsigned int n_refinements, ConvergenceTable &table) +{ + using Number = double; + using VectorType = LinearAlgebra::distributed::Vector; + + Triangulation tria; + + GridGenerator::hyper_cube(tria); + tria.refine_global(n_refinements); + + const MappingQ1 mapping; + const FE_Q fe_q(degree); + const FESystem fe(fe_q, n_components); + const QGauss quadrature(degree + 1); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + DoFTools::make_zero_boundary_constraints(dof_handler, constraints); + constraints.close(); + + for (unsigned int i = 0; i < 3; ++i) + { + LaplaceOperator + laplace_operator(i); + + laplace_operator.reinit(mapping, + dof_handler, + constraints, + quadrature.get_tensor_basis()[0]); + + VectorType src, dst; + + laplace_operator.initialize_dof_vector(src); + laplace_operator.initialize_dof_vector(dst); + + { + LinearAlgebra::distributed::Vector src_host( + src.get_partitioner()); + + VectorTools::create_right_hand_side( + mapping, + dof_handler, + quadrature, + AnalyticalFunction(n_components), + src_host, + constraints); + + LinearAlgebra::ReadWriteVector rw_vector( + src.get_partitioner()->locally_owned_range()); + rw_vector.import(src_host, VectorOperation::insert); + src.import(rw_vector, VectorOperation::insert); + + dst = 0.0; + } + + laplace_operator.vmult(dst, src); + + table.add_value("fe_degree", degree); + table.add_value("n_refinements", n_refinements); + table.add_value("n_components", n_components); + table.add_value("n_dofs", dof_handler.n_dofs()); + + if (std::is_same_v) + table.add_value("version", "host"); + else + table.add_value("version", "default"); + + table.add_value("norm", dst.l2_norm()); + table.set_scientific("norm", true); + } +} + +int +main() +{ + initlog(); + Kokkos::initialize(); + + const unsigned int dim = 2; + const unsigned int fe_degree = 3; + unsigned int n_refinements = 3; + + ConvergenceTable table; + + run(n_refinements, table); + run(n_refinements, table); + run(n_refinements, table); + run(n_refinements, table); + run(n_refinements, table); + run(n_refinements, table); + + table.write_text(deallog.get_file_stream()); + + Kokkos::finalize(); +} diff --git a/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.output b/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.output new file mode 100644 index 0000000000..d06c47ddd1 --- /dev/null +++ b/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.output @@ -0,0 +1,20 @@ + +fe_degree n_refinements n_components n_dofs version norm +3 3 1 625 host 1.1706e-04 +3 3 1 625 host 1.2457e-01 +3 3 1 625 host 1.2465e-01 +3 3 1 625 default 1.1706e-04 +3 3 1 625 default 1.2457e-01 +3 3 1 625 default 1.2465e-01 +3 3 2 1250 host 2.6174e-04 +3 3 2 1250 host 2.7854e-01 +3 3 2 1250 host 2.7872e-01 +3 3 2 1250 default 2.6174e-04 +3 3 2 1250 default 2.7854e-01 +3 3 2 1250 default 2.7872e-01 +3 3 3 1875 host 4.3798e-04 +3 3 3 1875 host 4.6609e-01 +3 3 3 1875 host 4.6639e-01 +3 3 3 1875 default 4.3798e-04 +3 3 3 1875 default 4.6609e-01 +3 3 3 1875 default 4.6639e-01 -- 2.39.5