Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3020-gf6ff73b199 2025-04-05 03:00:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
portable_matrix_free.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
16#ifndef dealii_portable_matrix_free_h
17#define dealii_portable_matrix_free_h
18
19#include <deal.II/base/config.h>
20
26#include <deal.II/base/tensor.h>
28
30
32#include <deal.II/fe/mapping.h>
33
35
38
40
41#include <Kokkos_Array.hpp>
42#include <Kokkos_Core.hpp>
43
44
45
47
48// Forward declaration
49namespace internal
50{
51 namespace MatrixFreeFunctions
52 {
53 enum class ConstraintKinds : std::uint16_t;
54 }
55} // namespace internal
56
57namespace Portable
58{
59 // forward declaration
60#ifndef DOXYGEN
61 namespace internal
62 {
63 template <int dim, typename Number>
64 class ReinitHelper;
65 }
66 template <int dim, typename Number>
67 struct SharedData;
68#endif
69
76 inline constexpr unsigned int n_max_dof_handlers = 5;
77
85 template <typename Number>
87 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>;
88
95 template <typename Number>
97 {
98 public:
103 : n_blocks(0)
104 {}
105
109 DeviceBlockVector(const DeviceBlockVector &other) = default;
110
116 : n_blocks(1)
117 , blocks{src}
118 {}
119
124 template <typename MemorySpace>
127 : n_blocks(src.n_blocks())
128 {
130 ExcMessage("Portable::MatrixFree is configured with " +
132 " but you are passing a BlockVector with " +
133 Utilities::to_string(src.n_blocks()) + " blocks."));
134
135 for (int b = 0; b < src.n_blocks(); ++b)
136 blocks[b] = DeviceVector<Number>(src.block(b).get_values(),
137 src.block(b).locally_owned_size());
138 }
139
145 block(unsigned int index)
146 {
148 return blocks[index];
149 }
150
155 block(unsigned int index) const
156 {
158 return blocks[index];
159 }
160
161 private:
165 unsigned int n_blocks;
166
170 Kokkos::Array<DeviceVector<Number>, n_max_dof_handlers> blocks;
171 };
172
197 template <int dim, typename Number = double>
199 {
200 public:
205
210 {
217 const bool use_coloring = false,
218 const bool overlap_communication_computation = false)
222 {
223#ifndef DEAL_II_MPI_WITH_DEVICE_SUPPORT
227 "Overlapping communication and computation requires device-aware MPI."));
228#endif
233 "Overlapping communication and coloring are incompatible options. Only one of them can be enabled."));
234 }
245
252
258 };
259
266 {
270 Kokkos::View<point_type **, MemorySpace::Default::kokkos_space> q_points;
271
276 Kokkos::View<types::global_dof_index **,
279
283 Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>
285
289 Kokkos::View<Number **, MemorySpace::Default::kokkos_space> JxW;
290
297
301 Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_values;
302
306 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
308
312 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
314
318 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
320
324 unsigned int n_cells;
325
329 unsigned int n_components;
330
334 unsigned int padding_length;
335
339 unsigned int row_start;
340
346
353
357 unsigned int scratch_pad_size;
358 };
359
360
368 struct Data
369 {
370 using TeamHandle = Kokkos::TeamPolicy<
371 MemorySpace::Default::kokkos_space::execution_space>::member_type;
372
377
378 const unsigned int n_dofhandler;
379 const int cell_index;
382
387 DEAL_II_HOST_DEVICE unsigned int
388 local_q_point_id(const unsigned int cell,
389 const unsigned int n_q_points,
390 const unsigned int q_point) const
391 {
393 cell) *
394 n_q_points +
395 q_point;
396 }
397
398
404 get_quadrature_point(const unsigned int cell,
405 const unsigned int q_point) const
406 {
407 return precomputed_data->q_points(q_point, cell);
408 }
409 };
410
415
419 unsigned int
421
432 template <typename IteratorFiltersType>
433 void
434 reinit(const Mapping<dim> &mapping,
436 const AffineConstraints<Number> &constraints,
437 const Quadrature<1> &quad,
438 const IteratorFiltersType &iterator_filter,
439 const AdditionalData &additional_data = AdditionalData());
440
444 void
445 reinit(const Mapping<dim> &mapping,
447 const AffineConstraints<Number> &constraints,
448 const Quadrature<1> &quad,
449 const AdditionalData &additional_data = AdditionalData());
450
454 void
456 const AffineConstraints<Number> &constraints,
457 const Quadrature<1> &quad,
458 const AdditionalData &additional_data = AdditionalData());
459
464 get_data(unsigned int color) const;
465
466 // clang-format off
481 // clang-format on
482 template <typename Functor, typename VectorType>
483 void
484 cell_loop(const Functor &func,
485 const VectorType &src,
486 VectorType &dst) const;
487
501 template <typename Functor>
502 void
503 evaluate_coefficients(Functor func) const;
504
509 template <typename VectorType>
510 void
511 copy_constrained_values(const VectorType &src, VectorType &dst) const;
512
518 template <typename VectorType>
519 void
520 set_constrained_values(const Number val, VectorType &dst) const;
521
527 template <typename MemorySpaceType>
528 void
531
535 const std::vector<std::vector<CellFilter>> &
537
548 const std::shared_ptr<const Utilities::MPI::Partitioner> &
550
554 const DoFHandler<dim> &
556
560 std::size_t
562
563 private:
567 template <typename IteratorFiltersType>
568 void
571 const AffineConstraints<Number> &constraints,
572 const Quadrature<1> &quad,
573 const IteratorFiltersType &iterator_filter,
574 const std::shared_ptr<const MPI_Comm> &comm,
575 const AdditionalData additional_data);
576
581 template <typename Functor, typename VectorType>
582 void
583 serial_cell_loop(const Functor &func,
584 const VectorType &src,
585 VectorType &dst) const;
586
591 template <typename Functor>
592 void
594 const Functor &func,
596 &src,
598 const;
599
603 template <typename Functor>
604 void
606 const Functor &func,
610 &dst) const;
611
612
616 int my_id;
617
624
630
635
642
646 unsigned int scratch_pad_size;
647
651 unsigned int fe_degree;
652
656 unsigned int n_components;
657
662
666 unsigned int dofs_per_cell;
667
671 unsigned int n_constrained_dofs;
672
676 unsigned int q_points_per_cell;
677
681 unsigned int n_colors;
682
686 std::vector<unsigned int> n_cells;
687
692 std::vector<Kokkos::View<point_type **, MemorySpace::Default::kokkos_space>>
694
699 std::vector<Kokkos::View<types::global_dof_index **,
702
707 std::vector<
708 Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>>
710
715 std::vector<Kokkos::View<Number **, MemorySpace::Default::kokkos_space>>
717
721 Kokkos::View<types::global_dof_index *, MemorySpace::Default::kokkos_space>
723
727 std::vector<
731
735 Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_values;
736
740 Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_gradients;
741
745 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
747
751 Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
753
758 std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
759
760
765 unsigned int padding_length;
766
770 std::vector<unsigned int> row_start;
771
776
780 std::vector<std::vector<CellFilter>> graph;
781
782 friend class internal::ReinitHelper<dim, Number>;
783 };
784
785
786
787 template <int dim, typename Number>
789 {
790 using SharedViewValues = Kokkos::View<
791 Number **,
792 MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
793 Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
794 using SharedViewGradients = Kokkos::View<
795 Number ***,
796 MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
797 Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
798 using SharedViewScratchPad = Kokkos::View<
799 Number *,
800 MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
801 Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
802
811
816
821
826 };
827
828
829
834 template <int dim, typename Number>
835 struct DataHost
836 {
840 typename Kokkos::View<Point<dim, Number> **,
843
848 typename Kokkos::View<types::global_dof_index **,
851
855 typename Kokkos::View<Number **[dim][dim],
858
862 typename Kokkos::View<Number **,
864
868 unsigned int n_cells;
869
873 unsigned int padding_length;
874
878 unsigned int row_start;
879
883 typename Kokkos::View<
886
892
893
894
898 unsigned int
899 local_q_point_id(const unsigned int cell,
900 const unsigned int n_q_points,
901 const unsigned int q_point) const
902 {
903 return (row_start / padding_length + cell) * n_q_points + q_point;
904 }
905
906
907
912 get_quadrature_point(const unsigned int cell,
913 const unsigned int q_point) const
914 {
915 return q_points(q_point, cell);
916 }
917 };
918
919
920
927 template <int dim, typename Number>
928 DataHost<dim, Number>
930 const typename ::Portable::MatrixFree<dim, Number>::PrecomputedData
931 &data,
932 const UpdateFlags &update_flags)
933 {
934 DataHost<dim, Number> data_host;
935
936 data_host.n_cells = data.n_cells;
937 data_host.padding_length = data.padding_length;
938 data_host.row_start = data.row_start;
939 data_host.use_coloring = data.use_coloring;
940
941 if (update_flags & update_quadrature_points)
942 {
943 data_host.q_points = Kokkos::create_mirror(data.q_points);
944 Kokkos::deep_copy(data_host.q_points, data.q_points);
945 }
946
947 data_host.local_to_global = Kokkos::create_mirror(data.local_to_global);
948 Kokkos::deep_copy(data_host.local_to_global, data.local_to_global);
949
950 if (update_flags & update_gradients)
951 {
952 data_host.inv_jacobian = Kokkos::create_mirror(data.inv_jacobian);
953 Kokkos::deep_copy(data_host.inv_jacobian, data.inv_jacobian);
954 }
955
956 if (update_flags & update_JxW_values)
957 {
958 data_host.JxW = Kokkos::create_mirror(data.JxW);
959 Kokkos::deep_copy(data_host.JxW, data.JxW);
960 }
961
962 data_host.constraint_mask = Kokkos::create_mirror(data.constraint_mask);
963 Kokkos::deep_copy(data_host.constraint_mask, data.constraint_mask);
964
965 return data_host;
966 }
967
968
969 /*----------------------- Inline functions ---------------------------------*/
970
971#ifndef DOXYGEN
972
973 template <int dim, typename Number>
974 inline const std::vector<std::vector<
977 {
978 return graph;
979 }
980
981
982
983 template <int dim, typename Number>
984 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
986 {
987 return partitioner;
988 }
989
990
991
992 template <int dim, typename Number>
993 inline const DoFHandler<dim> &
995 {
996 Assert(dof_handler != nullptr, ExcNotInitialized());
997
998 return *dof_handler;
999 }
1000
1001#endif
1002
1003} // namespace Portable
1004
1006
1007#endif
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
Abstract base class for mapping classes.
Definition mapping.h:320
DeviceBlockVector(const DeviceBlockVector &other)=default
Kokkos::Array< DeviceVector< Number >, n_max_dof_handlers > blocks
const DeviceVector< Number > & block(unsigned int index) const
DeviceBlockVector(const DeviceVector< Number > &src)
DeviceVector< Number > & block(unsigned int index)
DeviceBlockVector(const LinearAlgebra::distributed::BlockVector< Number, MemorySpace > &src)
const std::vector< std::vector< CellFilter > > & get_colored_graph() const
std::vector< std::vector< CellFilter > > graph
void reinit(const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const AdditionalData &additional_data=AdditionalData())
void internal_reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const IteratorFiltersType &iterator_filter, const std::shared_ptr< const MPI_Comm > &comm, const AdditionalData additional_data)
std::vector< Kokkos::View<::internal::MatrixFreeFunctions::ConstraintKinds *, MemorySpace::Default::kokkos_space > > constraint_mask
unsigned int get_padding_length() const
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights
const DoFHandler< dim > * dof_handler
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_gradients
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number, MemorySpaceType > &vec) const
void reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const IteratorFiltersType &iterator_filter, const AdditionalData &additional_data=AdditionalData())
std::vector< Kokkos::View< Number **[dim][dim], MemorySpace::Default::kokkos_space > > inv_jacobian
PrecomputedData get_data(unsigned int color) const
Kokkos::View< types::global_dof_index *, MemorySpace::Default::kokkos_space > constrained_dofs
std::vector< Kokkos::View< point_type **, MemorySpace::Default::kokkos_space > > q_points
types::global_dof_index n_dofs
void distributed_cell_loop(const Functor &func, const LinearAlgebra::distributed::Vector< Number, MemorySpace::Default > &src, LinearAlgebra::distributed::Vector< Number, MemorySpace::Default > &dst) const
std::vector< unsigned int > row_start
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > co_shape_gradients
std::vector< Kokkos::View< types::global_dof_index **, MemorySpace::Default::kokkos_space > > local_to_global
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_values
std::size_t memory_consumption() const
void set_constrained_values(const Number val, VectorType &dst) const
void distributed_cell_loop(const Functor &func, const LinearAlgebra::distributed::BlockVector< Number, MemorySpace::Default > &src, LinearAlgebra::distributed::BlockVector< Number, MemorySpace::Default > &dst) const
void copy_constrained_values(const VectorType &src, VectorType &dst) const
std::vector< Kokkos::View< Number **, MemorySpace::Default::kokkos_space > > JxW
const DoFHandler< dim > & get_dof_handler() const
std::shared_ptr< const Utilities::MPI::Partitioner > partitioner
void serial_cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
void evaluate_coefficients(Functor func) const
DataHost< dim, Number > copy_mf_data_to_host(const typename::Portable::MatrixFree< dim, Number >::PrecomputedData &data, const UpdateFlags &update_flags)
::internal::MatrixFreeFunctions::ElementType element_type
void reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const AdditionalData &additional_data=AdditionalData())
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner() const
std::vector< unsigned int > n_cells
void cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_HOST_DEVICE
Definition config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
std::vector< index_type > data
Definition mpi.cc:746
const MPI_Comm comm
Definition mpi.cc:924
constexpr unsigned int n_max_dof_handlers
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > DeviceVector
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:475
::Kokkos::DefaultExecutionSpace::memory_space kokkos_space
Kokkos::View< types::global_dof_index **, MemorySpace::Default::kokkos_space >::HostMirror local_to_global
Kokkos::View< Number **, MemorySpace::Default::kokkos_space >::HostMirror JxW
Point< dim, Number > get_quadrature_point(const unsigned int cell, const unsigned int q_point) const
Kokkos::View<::internal::MatrixFreeFunctions::ConstraintKinds *, MemorySpace::Default::kokkos_space >::HostMirror constraint_mask
Kokkos::View< Number **[dim][dim], MemorySpace::Default::kokkos_space >::HostMirror inv_jacobian
Kokkos::View< Point< dim, Number > **, MemorySpace::Default::kokkos_space >::HostMirror q_points
unsigned int local_q_point_id(const unsigned int cell, const unsigned int n_q_points, const unsigned int q_point) const
AdditionalData(const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values|update_quadrature_points, const bool use_coloring=false, const bool overlap_communication_computation=false)
unsigned int local_q_point_id(const unsigned int cell, const unsigned int n_q_points, const unsigned int q_point) const
SharedData< dim, Number > * shared_data
Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type TeamHandle
const PrecomputedData * precomputed_data
Portable::MatrixFree< dim, Number >::point_type & get_quadrature_point(const unsigned int cell, const unsigned int q_point) const
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_gradients
Kokkos::View< point_type **, MemorySpace::Default::kokkos_space > q_points
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_values
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > co_shape_gradients
Kokkos::View< Number **[dim][dim], MemorySpace::Default::kokkos_space > inv_jacobian
::internal::MatrixFreeFunctions::ElementType element_type
Kokkos::View< Number **, MemorySpace::Default::kokkos_space > JxW
Kokkos::View<::internal::MatrixFreeFunctions::ConstraintKinds *, MemorySpace::Default::kokkos_space > constraint_mask
Kokkos::View< types::global_dof_index **, MemorySpace::Default::kokkos_space > local_to_global
Kokkos::View< Number ***, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > SharedViewGradients
SharedViewScratchPad scratch_pad
Kokkos::View< Number *, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > SharedViewScratchPad
SharedData(const SharedViewValues &values, const SharedViewGradients &gradients, const SharedViewScratchPad &scratch_pad)
Kokkos::View< Number **, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > SharedViewValues
SharedViewGradients gradients