Reference documentation for deal.II version GIT 01a9543f1b 2023-12-05 20:40:02+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\}}\)
cuda_matrix_free.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #ifndef dealii_cuda_matrix_free_h
18 #define dealii_cuda_matrix_free_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/cuda_size.h>
24 #include <deal.II/base/mpi_stub.h>
27 #include <deal.II/base/tensor.h>
28 #include <deal.II/base/utilities.h>
29 
31 
33 #include <deal.II/fe/mapping.h>
34 
36 
40 
41 #include <Kokkos_Core.hpp>
42 
43 
44 
46 
47 // Forward declaration
48 namespace internal
49 {
50  namespace MatrixFreeFunctions
51  {
52  enum class ConstraintKinds : std::uint16_t;
53  }
54 } // namespace internal
55 
56 namespace CUDAWrappers
57 {
58  // forward declaration
59 #ifndef DOXYGEN
60  namespace internal
61  {
62  template <int dim, typename Number>
63  class ReinitHelper;
64  }
65 #endif
66 
93  template <int dim, typename Number = double>
94  class MatrixFree : public Subscriptor
95  {
96  public:
99  using CellFilter =
101 
106  {
110  AdditionalData(const UpdateFlags mapping_update_flags =
113  const bool use_coloring = false,
114  const bool overlap_communication_computation = false)
115  : mapping_update_flags(mapping_update_flags)
116  , use_coloring(use_coloring)
117  , overlap_communication_computation(overlap_communication_computation)
118  {
119 #ifndef DEAL_II_MPI_WITH_DEVICE_SUPPORT
120  AssertThrow(
121  overlap_communication_computation == false,
122  ExcMessage(
123  "Overlapping communication and computation requires CUDA-aware MPI."));
124 #endif
125  if (overlap_communication_computation == true)
126  AssertThrow(
127  use_coloring == false || overlap_communication_computation == false,
128  ExcMessage(
129  "Overlapping communication and coloring are incompatible options. Only one of them can be enabled."));
130  }
141 
148 
154  };
155 
160  struct Data
161  {
165  Kokkos::View<point_type **, MemorySpace::Default::kokkos_space> q_points;
166 
174 
178  Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>
180 
184  Kokkos::View<Number **, MemorySpace::Default::kokkos_space> JxW;
185 
192 
196  Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_values;
197 
201  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
203 
207  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
209 
213  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
215 
219  unsigned int n_cells;
220 
224  unsigned int padding_length;
225 
229  unsigned int row_start;
230 
236 
241  DEAL_II_HOST_DEVICE unsigned int
242  local_q_point_id(const unsigned int cell,
243  const unsigned int n_q_points,
244  const unsigned int q_point) const
245  {
246  return (row_start / padding_length + cell) * n_q_points + q_point;
247  }
248 
249 
255  get_quadrature_point(const unsigned int cell,
256  const unsigned int q_point) const
257  {
258  return q_points(cell, q_point);
259  }
260  };
261 
266 
270  unsigned int
272 
283  template <typename IteratorFiltersType>
284  void
285  reinit(const Mapping<dim> &mapping,
286  const DoFHandler<dim> &dof_handler,
287  const AffineConstraints<Number> &constraints,
288  const Quadrature<1> &quad,
289  const IteratorFiltersType &iterator_filter,
290  const AdditionalData &additional_data = AdditionalData());
291 
295  void
296  reinit(const Mapping<dim> &mapping,
297  const DoFHandler<dim> &dof_handler,
298  const AffineConstraints<Number> &constraints,
299  const Quadrature<1> &quad,
300  const AdditionalData &additional_data = AdditionalData());
301 
305  void
306  reinit(const DoFHandler<dim> &dof_handler,
307  const AffineConstraints<Number> &constraints,
308  const Quadrature<1> &quad,
309  const AdditionalData &additional_data = AdditionalData());
310 
314  Data
315  get_data(unsigned int color) const;
316 
317  // clang-format off
335  // clang-format on
336  template <typename Functor, typename VectorType>
337  void
338  cell_loop(const Functor &func,
339  const VectorType &src,
340  VectorType &dst) const;
341 
357  template <typename Functor>
358  void
359  evaluate_coefficients(Functor func) const;
360 
365  template <typename VectorType>
366  void
367  copy_constrained_values(const VectorType &src, VectorType &dst) const;
368 
374  template <typename VectorType>
375  void
376  set_constrained_values(const Number val, VectorType &dst) const;
377 
378 #ifdef DEAL_II_WITH_CUDA
383  void
386 #endif
387 
393  void
396  const;
397 
401  const std::vector<std::vector<CellFilter>> &
403 
414  const std::shared_ptr<const Utilities::MPI::Partitioner> &
416 
420  const DoFHandler<dim> &
422 
426  std::size_t
428 
429  private:
433  template <typename IteratorFiltersType>
434  void
436  const DoFHandler<dim> &dof_handler,
437  const AffineConstraints<Number> &constraints,
438  const Quadrature<1> &quad,
439  const IteratorFiltersType &iterator_filter,
440  const std::shared_ptr<const MPI_Comm> &comm,
441  const AdditionalData additional_data);
442 
447  template <typename Functor, typename VectorType>
448  void
449  serial_cell_loop(const Functor &func,
450  const VectorType &src,
451  VectorType &dst) const;
452 
457  template <typename Functor>
458  void
460  const Functor &func,
462  &src,
464  const;
465 
466 #ifdef DEAL_II_WITH_CUDA
472  template <typename Functor>
473  void
475  const Functor &func,
478 #endif
479 
483  int my_id;
484 
491 
497 
502 
506  unsigned int fe_degree;
507 
511  unsigned int dofs_per_cell;
512 
516  unsigned int n_constrained_dofs;
517 
521  unsigned int q_points_per_cell;
522 
526  unsigned int n_colors;
527 
531  std::vector<unsigned int> n_cells;
532 
537  std::vector<Kokkos::View<point_type **, MemorySpace::Default::kokkos_space>>
539 
544  std::vector<Kokkos::View<types::global_dof_index **,
547 
552  std::vector<
553  Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>>
555 
560  std::vector<Kokkos::View<Number **, MemorySpace::Default::kokkos_space>>
562 
566  Kokkos::View<types::global_dof_index *, MemorySpace::Default::kokkos_space>
568 
572  std::vector<
576 
580  Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_values;
581 
585  Kokkos::View<Number *, MemorySpace::Default::kokkos_space> shape_gradients;
586 
590  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
592 
596  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
598 
603  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
604 
605 
610  unsigned int padding_length;
611 
615  std::vector<unsigned int> row_start;
616 
621 
625  std::vector<std::vector<CellFilter>> graph;
626 
627  friend class internal::ReinitHelper<dim, Number>;
628  };
629 
630 
631 
632  template <int dim, typename Number>
633  struct SharedData
634  {
635  using TeamHandle = Kokkos::TeamPolicy<
636  MemorySpace::Default::kokkos_space::execution_space>::member_type;
637 
639  Number *,
640  MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
641  Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
643  Number *[dim],
644  MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
645  Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
646 
648  SharedData(const TeamHandle &team_member,
649  const SharedView1D &values,
650  const SharedView2D &gradients)
651  : team_member(team_member)
652  , values(values)
654  {}
655 
660 
665 
670  };
671 
672 
673 
678  template <int dim, typename Number>
679  struct DataHost
680  {
684  typename Kokkos::View<Point<dim, Number> **,
687 
695 
699  typename Kokkos::View<Number **[dim][dim],
702 
706  typename Kokkos::View<Number **,
708 
712  unsigned int n_cells;
713 
717  unsigned int padding_length;
718 
722  unsigned int row_start;
723 
727  typename Kokkos::View<
730 
736 
737 
738 
742  unsigned int
743  local_q_point_id(const unsigned int cell,
744  const unsigned int n_q_points,
745  const unsigned int q_point) const
746  {
747  return (row_start / padding_length + cell) * n_q_points + q_point;
748  }
749 
750 
751 
756  get_quadrature_point(const unsigned int cell,
757  const unsigned int q_point) const
758  {
759  return q_points(cell, q_point);
760  }
761  };
762 
763 
764 
771  template <int dim, typename Number>
772  DataHost<dim, Number>
774  const typename ::CUDAWrappers::MatrixFree<dim, Number>::Data &data,
775  const UpdateFlags &update_flags)
776  {
777  DataHost<dim, Number> data_host;
778 
779  data_host.n_cells = data.n_cells;
780  data_host.padding_length = data.padding_length;
781  data_host.row_start = data.row_start;
782  data_host.use_coloring = data.use_coloring;
783 
784  if (update_flags & update_quadrature_points)
785  {
786  data_host.q_points = Kokkos::create_mirror(data.q_points);
787  Kokkos::deep_copy(data_host.q_points, data.q_points);
788  }
789 
790  data_host.local_to_global = Kokkos::create_mirror(data.local_to_global);
791  Kokkos::deep_copy(data_host.local_to_global, data.local_to_global);
792 
793  if (update_flags & update_gradients)
794  {
795  data_host.inv_jacobian = Kokkos::create_mirror(data.inv_jacobian);
796  Kokkos::deep_copy(data_host.inv_jacobian, data.inv_jacobian);
797  }
798 
799  if (update_flags & update_JxW_values)
800  {
801  data_host.JxW = Kokkos::create_mirror(data.JxW);
802  Kokkos::deep_copy(data_host.JxW, data.JxW);
803  }
804 
805  data_host.constraint_mask = Kokkos::create_mirror(data.constraint_mask);
806  Kokkos::deep_copy(data_host.constraint_mask, data.constraint_mask);
807 
808  return data_host;
809  }
810 
811 
812  /*----------------------- Inline functions ---------------------------------*/
813 
814 #ifndef DOXYGEN
815 
816  template <int dim, typename Number>
817  inline const std::vector<std::vector<
820  {
821  return graph;
822  }
823 
824 
825 
826  template <int dim, typename Number>
827  inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
829  {
830  return partitioner;
831  }
832 
833 
834 
835  template <int dim, typename Number>
836  inline const DoFHandler<dim> &
838  {
839  Assert(dof_handler != nullptr, ExcNotInitialized());
840 
841  return *dof_handler;
842  }
843 
844 #endif
845 
846 } // namespace CUDAWrappers
847 
849 
850 #endif
std::vector< Kokkos::View< point_type **, MemorySpace::Default::kokkos_space > > q_points
void evaluate_coefficients(Functor func) const
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights
const DoFHandler< dim > * dof_handler
Data get_data(unsigned int color) const
unsigned int get_padding_length() const
const std::vector< std::vector< CellFilter > > & get_colored_graph() const
void initialize_dof_vector(LinearAlgebra::CUDAWrappers::Vector< Number > &vec) const
void reinit(const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const AdditionalData &additional_data=AdditionalData())
std::vector< unsigned int > n_cells
const DoFHandler< dim > & get_dof_handler() const
std::vector< Kokkos::View< Number **[dim][dim], MemorySpace::Default::kokkos_space > > inv_jacobian
std::vector< Kokkos::View< Number **, MemorySpace::Default::kokkos_space > > JxW
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner() const
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > co_shape_gradients
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_gradients
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())
void serial_cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
std::size_t memory_consumption() const
void distributed_cell_loop(const Functor &func, const LinearAlgebra::distributed::Vector< Number, MemorySpace::Default > &src, LinearAlgebra::distributed::Vector< Number, MemorySpace::Default > &dst) const
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_values
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number, MemorySpace::Default > &vec) const
Kokkos::View< types::global_dof_index *, MemorySpace::Default::kokkos_space > constrained_dofs
std::vector< std::vector< CellFilter > > graph
std::vector< Kokkos::View< types::global_dof_index **, MemorySpace::Default::kokkos_space > > local_to_global
void copy_constrained_values(const VectorType &src, VectorType &dst) const
std::vector< Kokkos::View<::internal::MatrixFreeFunctions::ConstraintKinds *, MemorySpace::Default::kokkos_space > > constraint_mask
std::shared_ptr< const Utilities::MPI::Partitioner > partitioner
void reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const AdditionalData &additional_data=AdditionalData())
types::global_dof_index n_dofs
void distributed_cell_loop(const Functor &func, const LinearAlgebra::CUDAWrappers::Vector< Number > &src, LinearAlgebra::CUDAWrappers::Vector< Number > &dst) const
std::vector< unsigned int > row_start
void set_constrained_values(const Number val, VectorType &dst) const
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)
DataHost< dim, Number > copy_mf_data_to_host(const typename::CUDAWrappers::MatrixFree< dim, Number >::Data &data, const UpdateFlags &update_flags)
void cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
Abstract base class for mapping classes.
Definition: mapping.h:317
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner(const unsigned int dof_handler_index=0) const
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
Definition: tensor.h:516
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
unsigned int global_dof_index
Definition: types.h:82
#define DEAL_II_HOST_DEVICE
Definition: numbers.h:35
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< types::global_dof_index **, MemorySpace::Default::kokkos_space >::HostMirror local_to_global
unsigned int local_q_point_id(const unsigned int cell, const unsigned int n_q_points, const unsigned int q_point) const
Kokkos::View< Number **, MemorySpace::Default::kokkos_space >::HostMirror JxW
Kokkos::View< Number **[dim][dim], MemorySpace::Default::kokkos_space >::HostMirror inv_jacobian
Kokkos::View< Point< dim, Number > **, MemorySpace::Default::kokkos_space >::HostMirror q_points
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)
Kokkos::View< types::global_dof_index **, MemorySpace::Default::kokkos_space > local_to_global
Kokkos::View< Number **[dim][dim], MemorySpace::Default::kokkos_space > inv_jacobian
Kokkos::View<::internal::MatrixFreeFunctions::ConstraintKinds *, MemorySpace::Default::kokkos_space > constraint_mask
Kokkos::View< Number **, MemorySpace::Default::kokkos_space > JxW
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > co_shape_gradients
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights
unsigned int local_q_point_id(const unsigned int cell, const unsigned int n_q_points, const unsigned int q_point) const
CUDAWrappers::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< Number *, MemorySpace::Default::kokkos_space > shape_values
Kokkos::View< point_type **, MemorySpace::Default::kokkos_space > q_points
Kokkos::View< Number *, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > SharedView1D
SharedData(const TeamHandle &team_member, const SharedView1D &values, const SharedView2D &gradients)
Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type TeamHandle
Kokkos::View< Number *[dim], MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > SharedView2D
::Kokkos::DefaultExecutionSpace::memory_space kokkos_space
Definition: memory_space.h:46
const MPI_Comm comm