--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/portable_matrix_free.h>
+#include <deal.II/matrix_free/portable_tensor_product_kernels.h>
+
+
+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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ struct FEEvaluationImpl
+ {
+ using TeamHandle = Kokkos::TeamPolicy<
+ MemorySpace::Default::kokkos_space::execution_space>::member_type;
+ using SharedView = Kokkos::View<Number *,
+ MemorySpace::Default::kokkos_space::
+ execution_space::scratch_memory_space,
+ Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+
+ DEAL_II_HOST_DEVICE static void
+ evaluate(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const typename MatrixFree<dim, Number>::Data *data,
+ SharedData<dim, Number> *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<EvaluatorVariant::evaluate_general,
+ dim,
+ fe_degree + 1,
+ n_q_points_1d,
+ Number>
+ 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<false>(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<dim, Number>::Data *data,
+ SharedData<dim, Number> *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<EvaluatorVariant::evaluate_general,
+ dim,
+ fe_degree + 1,
+ n_q_points_1d,
+ Number>
+ 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<false>(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<false>(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 <int dim, int fe_degree, typename Number>
+ struct FEEvaluationImplCollocation
+ {
+ DEAL_II_HOST_DEVICE static void
+ evaluate(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const typename MatrixFree<dim, Number>::Data *data,
+ SharedData<dim, Number> *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<EvaluatorVariant::evaluate_general,
+ dim,
+ fe_degree + 1,
+ fe_degree + 1,
+ Number>
+ 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<dim, Number>::Data *data,
+ SharedData<dim, Number> *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<EvaluatorVariant::evaluate_general,
+ dim,
+ fe_degree + 1,
+ fe_degree + 1,
+ Number>
+ 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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ struct FEEvaluationImplTransformToCollocation
+ {
+ DEAL_II_HOST_DEVICE static void
+ evaluate(const unsigned int n_components,
+ const EvaluationFlags::EvaluationFlags evaluation_flag,
+ const typename MatrixFree<dim, Number>::Data *data,
+ SharedData<dim, Number> *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<EvaluatorVariant::evaluate_general,
+ dim,
+ fe_degree + 1,
+ n_q_points_1d,
+ Number>
+ 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<dim, Number>::Data *data,
+ const SharedData<dim, Number> *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<EvaluatorVariant::evaluate_general,
+ dim,
+ fe_degree + 1,
+ n_q_points_1d,
+ Number>
+ 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
#include <deal.II/base/utilities.h>
#include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/portable_evaluation_kernels.h>
#include <deal.II/matrix_free/portable_hanging_nodes_internal.h>
#include <deal.II/matrix_free/portable_matrix_free.h>
#include <deal.II/matrix_free/portable_matrix_free.templates.h>
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.
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)
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));
});
}
}
template <int dim,
int fe_degree,
int n_q_points_1d,
- int n_components_,
+ int n_components,
typename Number>
DEAL_II_HOST_DEVICE void
- FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::evaluate(
- const EvaluationFlags::EvaluationFlags evaluate_flag)
+ FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number>::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<dim, fe_degree, Number>::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<dim, fe_degree, n_q_points_1d, Number>::
+ evaluate(n_components, evaluation_flag, data, shared_data);
+ }
+ else
+ {
+ Kokkos::abort("The element type is not yet supported by the portable "
+ "matrix-free module.");
}
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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<false>(
- 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<dim, fe_degree, Number>::
+ 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<dim, fe_degree, n_q_points_1d, Number>::
+ integrate(n_components, integration_flag, data, shared_data);
+ }
+ else
+ {
+ Kokkos::abort("The element type is not yet supported by the portable "
+ "matrix-free module.");
}
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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);
n_components_,
Number>::value_type
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
- 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;
}
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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);
typename Number>
DEAL_II_HOST_DEVICE void
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
- 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];
}
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
get_gradient(int q_point) const
{
+ Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
gradient_type grad;
if constexpr (n_components_ == 1)
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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)
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/matrix_free/shape_info.h>
+
#include <Kokkos_Core.hpp>
*/
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.
*/
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.
*/
Number ***,
MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+ using SharedViewScratchPad = Kokkos::View<
+ Number *,
+ MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
+ Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
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)
{}
/**
* Memory for computed gradients in reference coordinate system.
*/
SharedViewGradients gradients;
+
+ /**
+ * Memory for temporary arrays required by evaluation and integration.
+ */
+ SharedViewScratchPad scratch_pad;
};
MemorySpace::Default::kokkos_space::execution_space::
scratch_memory_space,
Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+ using SharedViewScratchPad =
+ Kokkos::View<Number *,
+ MemorySpace::Default::kokkos_space::execution_space::
+ scratch_memory_space,
+ Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
ApplyKernel(Functor func,
const typename MatrixFree<dim, Number>::Data gpu_data,
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);
}
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<dim, Number> shared_data(team_member, values, gradients);
+ SharedData<dim, Number> shared_data(team_member,
+ values,
+ gradients,
+ scratch_pad);
func(team_member.league_rank(), &gpu_data, &shared_data, src, dst);
}
};
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;
}
+ 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 <int dim, typename Number>
template <typename IteratorFiltersType>
void
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<unsigned int>(
- std::ceil(dim * std::log2(fe_degree + 1.)));
+ padding_length = 1 << static_cast<unsigned int>(std::ceil(
+ dim * std::log2(static_cast<float>(n_q_points_1d))));
dofs_per_cell = fe.n_dofs_per_cell();
n_components = fe.n_components();
::dealii::internal::MatrixFreeFunctions::ShapeInfo<Number> 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<Number *, MemorySpace::Default::kokkos_space>(
+ /**
+ * Helper function that copies or adds the first N entries of src to
+ * dst, depending on the template argument "add".
+ */
+ template <bool add, typename ViewTypeIn, typename ViewTypeOut>
+ 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 <int n_q_points_1d,
+ template <int n_rows,
+ int n_columns,
+ int direction,
typename Number,
- int direction,
- bool dof_to_quad,
+ bool contract_over_rows,
bool add,
- bool in_place,
typename ViewTypeIn,
typename ViewTypeOut>
DEAL_II_HOST_DEVICE void
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 <int n_q_points_1d,
+ template <int n_rows,
+ int n_columns,
+ int direction,
typename Number,
- int direction,
- bool dof_to_quad,
+ bool contract_over_rows,
bool add,
- bool in_place,
typename ViewTypeIn,
typename ViewTypeOut>
DEAL_II_HOST_DEVICE void
{
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<Kokkos::Rank<2>, 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 <int n_q_points_1d,
+ template <int n_rows,
+ int n_columns,
+ int direction,
typename Number,
- int direction,
- bool dof_to_quad,
+ bool contract_over_rows,
bool add,
- bool in_place,
typename ViewTypeIn,
typename ViewTypeOut>
DEAL_II_HOST_DEVICE void
{
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<Kokkos::Rank<3>, 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 <int dim,
- int n_q_points_1d,
+ int n_rows,
+ int n_columns,
typename Number,
int direction,
- bool dof_to_quad,
+ bool contract_over_rows,
bool add,
- bool in_place,
typename ViewTypeIn,
typename ViewTypeOut>
DEAL_II_HOST_DEVICE void
{
#if KOKKOS_VERSION >= 40000
if constexpr (dim == 1)
- apply_1d<n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
+ apply_1d<n_rows, n_columns, direction, Number, contract_over_rows, add>(
team_member, shape_data, in, out);
if constexpr (dim == 2)
- apply_2d<n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
+ apply_2d<n_rows, n_columns, direction, Number, contract_over_rows, add>(
team_member, shape_data, in, out);
if constexpr (dim == 3)
- apply_3d<n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
+ apply_3d<n_rows, n_columns, direction, Number, contract_over_rows, add>(
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 <EvaluatorVariant variant,
int dim,
- int fe_degree,
- int n_q_points_1d,
+ int n_rows,
+ int n_columns,
typename Number>
struct EvaluatorTensorProduct
{};
* Internal evaluator for 1d-3d shape function using the tensor product form
* of the basis functions.
*/
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int dim, int n_rows, int n_columns, typename Number>
struct EvaluatorTensorProduct<evaluate_general,
dim,
- fe_degree,
- n_q_points_1d,
+ n_rows,
+ n_columns,
Number>
{
public:
using TeamHandle = Kokkos::TeamPolicy<
MemorySpace::Default::kokkos_space::execution_space>::member_type;
+ using SharedView = Kokkos::View<Number *,
+ MemorySpace::Default::kokkos_space::
+ execution_space::scratch_memory_space,
+ Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
+
DEAL_II_HOST_DEVICE
EvaluatorTensorProduct(
const TeamHandle &team_member,
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
shape_gradients,
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
- co_shape_gradients);
-
- /**
- * Evaluate the finite element function at the quadrature points.
- */
- template <typename ViewType>
- DEAL_II_HOST_DEVICE void
- evaluate_values(ViewType u);
-
- /**
- * Evaluate the gradients of the finite element function at the quadrature
- * points.
- */
- template <typename ViewTypeIn, typename ViewTypeOut>
- 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 <typename ViewType1, typename ViewType2>
- DEAL_II_HOST_DEVICE void
- evaluate_values_and_gradients(ViewType1 u, ViewType2 grad_u);
-
- /**
- * Helper function for integrate(). Integrate the finite element function.
- */
- template <typename ViewType>
- DEAL_II_HOST_DEVICE void
- integrate_values(ViewType u);
-
- /**
- * Helper function for integrate(). Integrate the gradients of the finite
- * element function.
- */
- template <bool add, typename ViewType1, typename ViewType2>
- 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 <typename ViewType1, typename ViewType2>
- 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
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.
*/
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
co_shape_gradients;
+
+ /**
+ * Temporary storage for in-place evaluations.
+ */
+ SharedView temp;
};
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int dim, int n_rows, int n_columns, typename Number>
DEAL_II_HOST_DEVICE
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::
+ EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
EvaluatorTensorProduct(
const TeamHandle &team_member,
Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_values,
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
shape_gradients,
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
- 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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int dim, int n_rows, int n_columns, typename Number>
template <int direction,
bool dof_to_quad,
bool add,
typename ViewTypeIn,
typename ViewTypeOut>
DEAL_II_HOST_DEVICE void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::values(const ViewTypeIn in,
- ViewTypeOut out) const
+ EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
+ values(const ViewTypeIn in, ViewTypeOut out) const
{
- apply<dim, n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
- team_member, shape_values, in, out);
+ if constexpr (in_place)
+ {
+ apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, false>(
+ team_member, shape_values, in, temp);
+
+ populate_view<add>(team_member, out, temp, out.extent(0));
+ }
+ else
+ apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, add>(
+ team_member, shape_values, in, out);
}
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int dim, int n_rows, int n_columns, typename Number>
template <int direction,
bool dof_to_quad,
bool add,
typename ViewTypeIn,
typename ViewTypeOut>
DEAL_II_HOST_DEVICE void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::gradients(const ViewTypeIn in,
- ViewTypeOut out) const
+ EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
+ gradients(const ViewTypeIn in, ViewTypeOut out) const
{
- apply<dim, n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
- team_member, shape_gradients, in, out);
+ if constexpr (in_place)
+ {
+ apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, false>(
+ team_member, shape_gradients, in, temp);
+
+ populate_view<add>(team_member, out, temp, out.extent(0));
+ }
+ else
+ apply<dim, n_rows, n_columns, Number, direction, dof_to_quad, add>(
+ team_member, shape_gradients, in, out);
}
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
+ template <int dim, int n_rows, int n_columns, typename Number>
template <int direction,
bool dof_to_quad,
bool add,
typename ViewTypeIn,
typename ViewTypeOut>
DEAL_II_HOST_DEVICE void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::co_gradients(const ViewTypeIn in,
- ViewTypeOut out) const
- {
- apply<dim, n_q_points_1d, Number, direction, dof_to_quad, add, in_place>(
- team_member, co_shape_gradients, in, out);
- }
-
-
-
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <typename ViewType>
- DEAL_II_HOST_DEVICE inline void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <typename ViewType>
- DEAL_II_HOST_DEVICE inline void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::integrate_values(ViewType u)
+ EvaluatorTensorProduct<evaluate_general, dim, n_rows, n_columns, Number>::
+ 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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <typename ViewTypeIn, typename ViewTypeOut>
- DEAL_II_HOST_DEVICE inline void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <typename ViewType1, typename ViewType2>
- DEAL_II_HOST_DEVICE inline void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <bool add, typename ViewType1, typename ViewType2>
- DEAL_II_HOST_DEVICE inline void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::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 <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <typename ViewType1, typename ViewType2>
- DEAL_II_HOST_DEVICE inline void
- EvaluatorTensorProduct<evaluate_general,
- dim,
- fe_degree,
- n_q_points_1d,
- Number>::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<dim,
+ n_columns,
+ n_columns,
+ Number,
+ direction,
+ dof_to_quad,
+ false>(team_member, co_shape_gradients, in, temp);
+
+ populate_view<add>(team_member, out, temp, out.extent(0));
}
else
- Kokkos::abort("dim must not exceed 3!");
+ apply<dim, n_columns, n_columns, Number, direction, dof_to_quad, add>(
+ team_member, co_shape_gradients, in, out);
}
} // namespace internal
} // namespace Portable
Kokkos::View<double *, MemorySpace::Default::kokkos_space> dst,
Kokkos::View<double *, MemorySpace::Default::kokkos_space> src)
{
+ Kokkos::View<
+ double *,
+ MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
+ Kokkos::MemoryTraits<Kokkos::Unmanaged>>
+ 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);
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");
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();
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
Kokkos::View<double *, MemorySpace::Default::kokkos_space> shape_gradients,
Kokkos::View<double *, MemorySpace::Default::kokkos_space> co_shape_gradients,
Kokkos::View<double *, MemorySpace::Default::kokkos_space> dst,
- Kokkos::View<double *, MemorySpace::Default::kokkos_space> src)
+ Kokkos::View<double *, MemorySpace::Default::kokkos_space> src,
+ Kokkos::View<double *, MemorySpace::Default::kokkos_space> tmp)
{
+ Kokkos::View<
+ double *,
+ MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
+ Kokkos::MemoryTraits<Kokkos::Unmanaged>>
+ 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<d0, dof_to_quad, false, false>(src, tmp);
team_member.team_barrier();
- evaluator.template values<1, dof_to_quad, add, false>(src, dst);
+ evaluator.template values<d1, dof_to_quad, add, false>(tmp, dst);
}
if (type == 1)
{
- evaluator.template gradients<0, dof_to_quad, false, false>(src, src);
+ evaluator.template gradients<d0, dof_to_quad, false, false>(src, tmp);
team_member.team_barrier();
- evaluator.template gradients<1, dof_to_quad, add, false>(src, dst);
+ evaluator.template gradients<d1, dof_to_quad, add, false>(tmp, dst);
}
}
constexpr int M_2d = M * M;
constexpr int N_2d = N * N;
+ constexpr int MN = M * N;
LinearAlgebra::ReadWriteVector<double> x_ref(N_2d), y_ref(M_2d);
Kokkos::View<double[N_2d], MemorySpace::Default::kokkos_space> x_dev(
Kokkos::view_alloc("x_dev", Kokkos::WithoutInitializing));
Kokkos::View<double[M_2d], MemorySpace::Default::kokkos_space> y_dev(
Kokkos::view_alloc("y_dev", Kokkos::WithoutInitializing));
+ Kokkos::View<double[MN], MemorySpace::Default::kokkos_space> 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);
shape_gradients,
co_shape_gradients,
y_dev,
- x_dev);
+ x_dev,
+ tmp_dev);
});
// Check the results on the host
shape_gradients,
co_shape_gradients,
x_dev,
- y_dev);
+ y_dev,
+ tmp_dev);
});
// Check the results on the host
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");
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();
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
// 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