--- /dev/null
+Improved: Portable::MatrixFree now supports FESystem(FE_Q, n_components)
+with arbitrary number of components.
+<br>
+(Peter Munch, Daniel Arndt, 2024/08/14)
/**
* 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<n_components_ == dim,
+ Tensor<2, dim, Number>,
+ Tensor<1, n_components_, Tensor<1, dim, Number>>>>;
/**
* An alias to kernel specific information.
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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<dim, fe_degree, false>(
- 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<dim, fe_degree, false, Number>(
+ shared_data->team_member,
+ data->constraint_weights,
+ data->constraint_mask(cell_id),
+ Kokkos::subview(shared_data->values, Kokkos::ALL, c));
+ }
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
distribute_local_to_global(Number *dst) const
{
- static_assert(n_components_ == 1, "This function only supports FE with one \
- components");
-
- internal::resolve_hanging_nodes<dim, fe_degree, true>(
- 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<dim, fe_degree, true, Number>(
+ 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));
});
}
}
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();
+ }
}
}
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<false>(
- 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<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();
+ }
}
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::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;
+ }
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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;
+ }
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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);
+ }
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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];
+ }
}
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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;
FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
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);
+ }
}
}
template <unsigned int fe_degree,
unsigned int direction,
bool transpose,
- typename Number>
+ typename Number,
+ typename ViewType>
DEAL_II_HOST_DEVICE inline void
interpolate_boundary_2d(
const Kokkos::TeamPolicy<
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
constraint_weights,
const dealii::internal::MatrixFreeFunctions::ConstraintKinds
- &constraint_mask,
- Kokkos::View<Number *,
- MemorySpace::Default::kokkos_space::execution_space::
- scratch_memory_space,
- Kokkos::MemoryTraits<Kokkos::Unmanaged>> 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);
template <unsigned int fe_degree,
unsigned int direction,
bool transpose,
- typename Number>
+ typename Number,
+ typename ViewType>
DEAL_II_HOST_DEVICE inline void
interpolate_boundary_3d(
const Kokkos::TeamPolicy<
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
constraint_weights,
const dealii::internal::MatrixFreeFunctions::ConstraintKinds
- constraint_mask,
- Kokkos::View<Number *,
- MemorySpace::Default::kokkos_space::execution_space::
- scratch_memory_space,
- Kokkos::MemoryTraits<Kokkos::Unmanaged>> 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);
* @cite ljungkvist2017matrix and in Section 3.4 of
* @cite kronbichler2019multigrid.
*/
- template <int dim, int fe_degree, bool transpose, typename Number>
+ template <int dim,
+ int fe_degree,
+ bool transpose,
+ typename Number,
+ typename ViewType>
DEAL_II_HOST_DEVICE void
resolve_hanging_nodes(
const Kokkos::TeamPolicy<
Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
constraint_weights,
const dealii::internal::MatrixFreeFunctions::ConstraintKinds
- constraint_mask,
- Kokkos::View<Number *,
- MemorySpace::Default::kokkos_space::execution_space::
- scratch_memory_space,
- Kokkos::MemoryTraits<Kokkos::Unmanaged>> values)
+ constraint_mask,
+ ViewType values)
{
if constexpr (dim == 2)
{
*/
unsigned int n_cells;
+ /**
+ * Number of components.
+ */
+ unsigned int n_components;
+
/**
* Length of the padding.
*/
*/
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.
*/
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<Kokkos::Unmanaged>>;
- using SharedView2D = Kokkos::View<
- Number *[dim],
+ using SharedViewGradients = 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 SharedView1D &values,
- const SharedView2D &gradients)
+ SharedData(const TeamHandle &team_member,
+ const SharedViewValues &values,
+ const SharedViewGradients &gradients)
: team_member(team_member)
, values(values)
, gradients(gradients)
/**
* Memory for dof and quad values.
*/
- SharedView1D values;
+ SharedViewValues values;
/**
* Memory for computed gradients in reference coordinate system.
*/
- SharedView2D gradients;
+ SharedViewGradients gradients;
};
const std::vector<unsigned int> &lexicographic_inv;
std::vector<types::global_dof_index> 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;
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)
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] =
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] =
{
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<Kokkos::Unmanaged>>;
- using SharedView2D =
- Kokkos::View<Number *[dim],
+ using SharedViewGradients =
+ Kokkos::View<Number ***,
MemorySpace::Default::kokkos_space::execution_space::
scratch_memory_space,
Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
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);
}
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<dim, Number> shared_data(team_member, values, gradients);
func(team_member.league_rank(), &gpu_data, &shared_data, src, dst);
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;
padding_length = 1 << static_cast<unsigned int>(
std::ceil(dim * std::log2(fe_degree + 1.)));
- dofs_per_cell = fe.n_dofs_per_cell();
- q_points_per_cell = Utilities::fixed_power<dim>(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<dim>(n_q_points_1d);
::dealii::internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info(quad,
fe);
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 <deal.II/base/convergence_table.h>
+
+#include <deal.II/distributed/fully_distributed_tria.h>
+#include <deal.II/distributed/repartitioning_policy_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/portable_fe_evaluation.h>
+#include <deal.II/matrix_free/portable_matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "Kokkos_Core.hpp"
+
+template <int dim,
+ int fe_degree,
+ int n_components,
+ typename Number,
+ typename MemorySpace>
+class LaplaceOperator;
+
+template <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperator<dim, fe_degree, n_components, Number, MemorySpace::Host>
+{
+public:
+ using VectorType =
+ LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>;
+
+ LaplaceOperator(const unsigned int version)
+ : version(version)
+ {}
+
+ void
+ reinit(const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof_handler,
+ const AffineConstraints<Number> &constraints,
+ const Quadrature<1> &quadrature)
+ {
+ typename MatrixFree<dim, Number>::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<dim, Number> &data,
+ VectorType &dst,
+ const VectorType &src,
+ const std::pair<unsigned int, unsigned int> &cell_range) const
+ {
+ FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, Number> 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<dim, Number> matrix_free;
+};
+
+
+
+template <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperatorQuad
+{
+public:
+ DEAL_II_HOST_DEVICE
+ LaplaceOperatorQuad(const unsigned int version)
+ : version(version)
+ {}
+
+ DEAL_II_HOST_DEVICE void
+ operator()(
+ Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, Number>
+ *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 <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperatorLocal
+{
+public:
+ LaplaceOperatorLocal(const unsigned int version)
+ : version(version)
+ {}
+
+ DEAL_II_HOST_DEVICE void
+ operator()(const unsigned int cell,
+ const typename Portable::MatrixFree<dim, Number>::Data *gpu_data,
+ Portable::SharedData<dim, Number> *shared_data,
+ const Number *src,
+ Number *dst) const
+ {
+ (void)cell; // TODO?
+
+ Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, Number>
+ 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<dim, fe_degree, n_components, Number>(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 <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperator<dim,
+ fe_degree,
+ n_components,
+ Number,
+ MemorySpace::Default>
+{
+public:
+ using VectorType =
+ LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>;
+
+ LaplaceOperator(const unsigned int version)
+ : version(version)
+ {}
+
+ void
+ reinit(const Mapping<dim> &mapping,
+ const DoFHandler<dim> &dof_handler,
+ const AffineConstraints<Number> &constraints,
+ const Quadrature<1> &quadrature)
+ {
+ typename Portable::MatrixFree<dim, Number>::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<dim, fe_degree, n_components, Number> 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<dim, Number> matrix_free;
+};
+
+
+
+template <int dim, typename T>
+class AnalyticalFunction : public Function<dim, T>
+{
+public:
+ AnalyticalFunction(const unsigned int n_components)
+ : Function<dim, T>(n_components)
+ {}
+
+ virtual T
+ value(const Point<dim, T> &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 <unsigned int dim,
+ const int degree,
+ int n_components,
+ typename MemorySpace>
+void
+run(const unsigned int n_refinements, ConvergenceTable &table)
+{
+ using Number = double;
+ using VectorType = LinearAlgebra::distributed::Vector<Number, MemorySpace>;
+
+ Triangulation<dim> tria;
+
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(n_refinements);
+
+ const MappingQ1<dim> mapping;
+ const FE_Q<dim> fe_q(degree);
+ const FESystem<dim> fe(fe_q, n_components);
+ const QGauss<dim> quadrature(degree + 1);
+
+ DoFHandler<dim> dof_handler(tria);
+ dof_handler.distribute_dofs(fe);
+
+ AffineConstraints<Number> constraints;
+ DoFTools::make_zero_boundary_constraints(dof_handler, constraints);
+ constraints.close();
+
+ for (unsigned int i = 0; i < 3; ++i)
+ {
+ LaplaceOperator<dim, degree, n_components, Number, MemorySpace>
+ 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<Number> src_host(
+ src.get_partitioner());
+
+ VectorTools::create_right_hand_side<dim, dim>(
+ mapping,
+ dof_handler,
+ quadrature,
+ AnalyticalFunction<dim, Number>(n_components),
+ src_host,
+ constraints);
+
+ LinearAlgebra::ReadWriteVector<Number> 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<MemorySpace, dealii::MemorySpace::Host>)
+ 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<dim, fe_degree, 1, MemorySpace::Host>(n_refinements, table);
+ run<dim, fe_degree, 1, MemorySpace::Default>(n_refinements, table);
+ run<dim, fe_degree, dim, MemorySpace::Host>(n_refinements, table);
+ run<dim, fe_degree, dim, MemorySpace::Default>(n_refinements, table);
+ run<dim, fe_degree, dim + 1, MemorySpace::Host>(n_refinements, table);
+ run<dim, fe_degree, dim + 1, MemorySpace::Default>(n_refinements, table);
+
+ table.write_text(deallog.get_file_stream());
+
+ Kokkos::finalize();
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff