Reference documentation for deal.II version GIT f8d87ffe95 2022-08-04 17:30: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 - 2022 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 #ifdef DEAL_II_COMPILER_CUDA_AWARE
23 
24 # include <deal.II/base/cuda_size.h>
25 # include <deal.II/base/mpi.h>
26 # include <deal.II/base/quadrature.h>
27 # include <deal.II/base/tensor.h>
28 
30 
32 # include <deal.II/fe/mapping.h>
33 # include <deal.II/fe/mapping_q1.h>
34 
36 
38 # include <deal.II/lac/cuda_vector.h>
40 
41 
42 
44 
45 // Forward declaration
46 namespace internal
47 {
48  namespace MatrixFreeFunctions
49  {
50  enum class ConstraintKinds : std::uint16_t;
51  }
52 } // namespace internal
53 
54 namespace CUDAWrappers
55 {
56  // forward declaration
57 # ifndef DOXYGEN
58  namespace internal
59  {
60  template <int dim, typename Number>
61  class ReinitHelper;
62  }
63 # endif
64 
91  template <int dim, typename Number = double>
92  class MatrixFree : public Subscriptor
93  {
94  public:
97  using CellFilter =
99 
106  {
108  parallel_over_elem
109  };
110 
115  {
120  const ParallelizationScheme parallelization_scheme = parallel_in_elem,
121  const UpdateFlags mapping_update_flags = update_gradients |
124  const bool use_coloring = false,
125  const bool overlap_communication_computation = false)
126  : parallelization_scheme(parallelization_scheme)
127  , mapping_update_flags(mapping_update_flags)
128  , use_coloring(use_coloring)
129  , overlap_communication_computation(overlap_communication_computation)
130  {
131 # ifndef DEAL_II_MPI_WITH_CUDA_SUPPORT
132  AssertThrow(
133  overlap_communication_computation == false,
134  ExcMessage(
135  "Overlapping communication and computation requires CUDA-aware MPI."));
136 # endif
137  if (overlap_communication_computation == true)
138  AssertThrow(
139  use_coloring == false || overlap_communication_computation == false,
140  ExcMessage(
141  "Overlapping communication and coloring are incompatible options. Only one of them can be enabled."));
142  }
158 
165 
171  };
172 
177  struct Data
178  {
183 
189 
193  Number *inv_jacobian;
194 
198  Number *JxW;
199 
203  unsigned int id;
204 
208  unsigned int n_cells;
209 
213  unsigned int padding_length;
214 
218  unsigned int row_start;
219 
224 
230  };
231 
236 
241 
245  unsigned int
247 
258  template <typename IteratorFiltersType>
259  void
260  reinit(const Mapping<dim> & mapping,
261  const DoFHandler<dim> & dof_handler,
262  const AffineConstraints<Number> &constraints,
263  const Quadrature<1> & quad,
264  const IteratorFiltersType & iterator_filter,
265  const AdditionalData & additional_data = AdditionalData());
266 
270  void
271  reinit(const Mapping<dim> & mapping,
272  const DoFHandler<dim> & dof_handler,
273  const AffineConstraints<Number> &constraints,
274  const Quadrature<1> & quad,
275  const AdditionalData & additional_data = AdditionalData());
276 
280  void
281  reinit(const DoFHandler<dim> & dof_handler,
282  const AffineConstraints<Number> &constraints,
283  const Quadrature<1> & quad,
284  const AdditionalData & additional_data = AdditionalData());
285 
289  Data
290  get_data(unsigned int color) const;
291 
292  // clang-format off
310  // clang-format on
311  template <typename Functor, typename VectorType>
312  void
313  cell_loop(const Functor & func,
314  const VectorType &src,
315  VectorType & dst) const;
316 
332  template <typename Functor>
333  void
334  evaluate_coefficients(Functor func) const;
335 
340  template <typename VectorType>
341  void
342  copy_constrained_values(const VectorType &src, VectorType &dst) const;
343 
349  template <typename VectorType>
350  void
351  set_constrained_values(const Number val, VectorType &dst) const;
352 
357  void
360 
366  void
369 
373  const std::vector<std::vector<CellFilter>> &
375 
386  const std::shared_ptr<const Utilities::MPI::Partitioner> &
388 
392  void
393  free();
394 
398  const DoFHandler<dim> &
400 
404  std::size_t
406 
407  private:
411  template <typename IteratorFiltersType>
412  void
413  internal_reinit(const Mapping<dim> & mapping,
414  const DoFHandler<dim> & dof_handler,
415  const AffineConstraints<Number> &constraints,
416  const Quadrature<1> & quad,
417  const IteratorFiltersType & iterator_filter,
418  std::shared_ptr<const MPI_Comm> comm,
419  const AdditionalData additional_data);
420 
425  template <typename Functor, typename VectorType>
426  void
427  serial_cell_loop(const Functor & func,
428  const VectorType &src,
429  VectorType & dst) const;
430 
435  template <typename Functor>
436  void
438  const Functor & func,
441 
447  template <typename Functor>
448  void
450  const Functor & func,
453 
458  template <typename VectorType>
459  void
460  serial_copy_constrained_values(const VectorType &src,
461  VectorType & dst) const;
462 
467  void
471 
478  void
482 
487  template <typename VectorType>
488  void
489  serial_set_constrained_values(const Number val, VectorType &dst) const;
490 
495  void
497  const Number val,
499 
506  void
508  const Number val,
510 
514  int my_id;
515 
521 
528 
534 
539 
543  unsigned int fe_degree;
544 
548  unsigned int dofs_per_cell;
549 
553  unsigned int n_constrained_dofs;
554 
558  unsigned int q_points_per_cell;
559 
563  unsigned int n_colors;
564 
568  std::vector<unsigned int> n_cells;
569 
574  std::vector<point_type *> q_points;
575 
580  std::vector<types::global_dof_index *> local_to_global;
581 
586  std::vector<Number *> inv_jacobian;
587 
592  std::vector<Number *> JxW;
593 
598 
602  std::vector<::internal::MatrixFreeFunctions::ConstraintKinds *>
604 
609  std::vector<dim3> grid_dim;
610 
615  std::vector<dim3> block_dim;
616 
621  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
622 
626  unsigned int cells_per_block;
627 
633 
639 
644  unsigned int padding_length;
645 
649  std::vector<unsigned int> row_start;
650 
655 
659  std::vector<std::vector<CellFilter>> graph;
660 
661  friend class internal::ReinitHelper<dim, Number>;
662  };
663 
664 
665 
666  // TODO find a better place to put these things
667 
671  template <int dim, typename Number>
672  struct SharedData
673  {
677  __device__
678  SharedData(Number *vd, Number *gq[dim])
679  : values(vd)
680  {
681  for (unsigned int d = 0; d < dim; ++d)
682  gradients[d] = gq[d];
683  }
684 
688  Number *values;
689 
695  Number *gradients[dim];
696  };
697 
698 
699 
700  // This function determines the number of cells per block, possibly at compile
701  // time (by virtue of being 'constexpr')
702  // TODO this function should be rewritten using meta-programming
703  __host__ __device__ constexpr unsigned int
704  cells_per_block_shmem(int dim, int fe_degree)
705  {
706  /* clang-format off */
707  // We are limiting the number of threads according to the
708  // following formulas:
709  // - in 2D: `threads = cells * (k+1)^d <= 4*CUDAWrappers::warp_size`
710  // - in 3D: `threads = cells * (k+1)^d <= 2*CUDAWrappers::warp_size`
711  return dim==2 ? (fe_degree==1 ? CUDAWrappers::warp_size : // 128
712  fe_degree==2 ? CUDAWrappers::warp_size/4 : // 72
713  fe_degree==3 ? CUDAWrappers::warp_size/8 : // 64
714  fe_degree==4 ? CUDAWrappers::warp_size/8 : // 100
715  1) :
716  dim==3 ? (fe_degree==1 ? CUDAWrappers::warp_size/4 : // 64
717  fe_degree==2 ? CUDAWrappers::warp_size/16 : // 54
718  1) : 1;
719  /* clang-format on */
720  }
721 
722 
723  /*----------------------- Helper functions ---------------------------------*/
729  template <int dim>
730  __device__ inline unsigned int
731  q_point_id_in_cell(const unsigned int n_q_points_1d)
732  {
733  return (
734  dim == 1 ? threadIdx.x % n_q_points_1d :
735  dim == 2 ? threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
736  threadIdx.x % n_q_points_1d +
737  n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
738  }
739 
740 
741 
748  template <int dim, typename Number>
749  __device__ inline unsigned int
751  const unsigned int cell,
752  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
753  const unsigned int n_q_points_1d,
754  const unsigned int n_q_points)
755  {
756  return (data->row_start / data->padding_length + cell) * n_q_points +
757  q_point_id_in_cell<dim>(n_q_points_1d);
758  }
759 
760 
761 
767  template <int dim, typename Number>
768  __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
770  const unsigned int cell,
771  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
772  const unsigned int n_q_points_1d)
773  {
774  return *(data->q_points + data->padding_length * cell +
775  q_point_id_in_cell<dim>(n_q_points_1d));
776  }
777 
782  template <int dim, typename Number>
783  struct DataHost
784  {
788  std::vector<Point<dim, Number>> q_points;
789 
794  std::vector<types::global_dof_index> local_to_global;
795 
799  std::vector<Number> inv_jacobian;
800 
804  std::vector<Number> JxW;
805 
809  unsigned int id;
810 
814  unsigned int n_cells;
815 
819  unsigned int padding_length;
820 
824  unsigned int row_start;
825 
829  std::vector<::internal::MatrixFreeFunctions::ConstraintKinds>
831 
837  };
838 
839 
840 
847  template <int dim, typename Number>
850  const typename ::CUDAWrappers::MatrixFree<dim, Number>::Data &data,
851  const UpdateFlags &update_flags)
852  {
853  DataHost<dim, Number> data_host;
854 
855  data_host.id = data.id;
856  data_host.n_cells = data.n_cells;
857  data_host.padding_length = data.padding_length;
858  data_host.row_start = data.row_start;
859  data_host.use_coloring = data.use_coloring;
860 
861  const unsigned int n_elements =
862  data_host.n_cells * data_host.padding_length;
863  if (update_flags & update_quadrature_points)
864  {
865  data_host.q_points.resize(n_elements);
866  Utilities::CUDA::copy_to_host(data.q_points, data_host.q_points);
867  }
868 
869  data_host.local_to_global.resize(n_elements);
870  Utilities::CUDA::copy_to_host(data.local_to_global,
871  data_host.local_to_global);
872 
873  if (update_flags & update_gradients)
874  {
875  data_host.inv_jacobian.resize(n_elements * dim * dim);
876  Utilities::CUDA::copy_to_host(data.inv_jacobian,
877  data_host.inv_jacobian);
878  }
879 
880  if (update_flags & update_JxW_values)
881  {
882  data_host.JxW.resize(n_elements);
883  Utilities::CUDA::copy_to_host(data.JxW, data_host.JxW);
884  }
885 
886  data_host.constraint_mask.resize(data_host.n_cells);
887  Utilities::CUDA::copy_to_host(data.constraint_mask,
888  data_host.constraint_mask);
889 
890  return data_host;
891  }
892 
893 
894 
900  template <int dim, typename Number>
901  inline unsigned int
902  local_q_point_id_host(const unsigned int cell,
903  const DataHost<dim, Number> &data,
904  const unsigned int n_q_points,
905  const unsigned int i)
906  {
907  return (data.row_start / data.padding_length + cell) * n_q_points + i;
908  }
909 
910 
911 
919  template <int dim, typename Number>
920  inline Point<dim, Number>
921  get_quadrature_point_host(const unsigned int cell,
922  const DataHost<dim, Number> &data,
923  const unsigned int i)
924  {
925  return data.q_points[data.padding_length * cell + i];
926  }
927 
928 
929  /*----------------------- Inline functions ---------------------------------*/
930 
931 # ifndef DOXYGEN
932 
933  template <int dim, typename Number>
934  inline const std::vector<std::vector<
937  {
938  return graph;
939  }
940 
941 
942 
943  template <int dim, typename Number>
944  inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
946  {
947  return partitioner;
948  }
949 
950 
951 
952  template <int dim, typename Number>
953  inline const DoFHandler<dim> &
955  {
956  Assert(dof_handler != nullptr, ExcNotInitialized());
957 
958  return *dof_handler;
959  }
960 
961 # endif
962 
963 } // namespace CUDAWrappers
964 
966 
967 #endif
968 #endif
void evaluate_coefficients(Functor func) const
void serial_copy_constrained_values(const VectorType &src, VectorType &dst) const
std::vector< Number * > JxW
const DoFHandler< dim > * dof_handler
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &vec) const
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 serial_set_constrained_values(const Number val, VectorType &dst) const
std::vector< dim3 > block_dim
void reinit(const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const AdditionalData &additional_data=AdditionalData())
void distributed_copy_constrained_values(const LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &src, LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &dst) const
std::vector< unsigned int > n_cells
const DoFHandler< dim > & get_dof_handler() const
ParallelizationScheme parallelization_scheme
void internal_reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const IteratorFiltersType &iterator_filter, std::shared_ptr< const MPI_Comm > comm, const AdditionalData additional_data)
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner() 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())
void serial_cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
std::size_t memory_consumption() const
std::vector<::internal::MatrixFreeFunctions::ConstraintKinds * > constraint_mask
std::vector< point_type * > q_points
CUDAWrappers::MatrixFree< dim, Number >::point_type & get_quadrature_point(const unsigned int cell, const typename CUDAWrappers::MatrixFree< dim, Number >::Data *data, const unsigned int n_q_points_1d)
void distributed_set_constrained_values(const Number val, LinearAlgebra::CUDAWrappers::Vector< Number > &dst) const
unsigned int local_q_point_id(const unsigned int cell, const typename CUDAWrappers::MatrixFree< dim, Number >::Data *data, const unsigned int n_q_points_1d, const unsigned int n_q_points)
void distributed_set_constrained_values(const Number val, LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &dst) const
Point< dim, Number > get_quadrature_point_host(const unsigned int cell, const DataHost< dim, Number > &data, const unsigned int i)
void distributed_cell_loop(const Functor &func, const LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &src, LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &dst) const
std::vector< dim3 > grid_dim
std::vector< types::global_dof_index * > local_to_global
std::vector< std::vector< CellFilter > > graph
void copy_constrained_values(const VectorType &src, VectorType &dst) const
std::shared_ptr< const Utilities::MPI::Partitioner > partitioner
unsigned int q_point_id_in_cell(const unsigned int n_q_points_1d)
void distributed_copy_constrained_values(const LinearAlgebra::CUDAWrappers::Vector< Number > &src, LinearAlgebra::CUDAWrappers::Vector< Number > &dst) const
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
DataHost< dim, Number > copy_mf_data_to_host(const typename::CUDAWrappers::MatrixFree< dim, Number >::Data &data, const UpdateFlags &update_flags)
types::global_dof_index * constrained_dofs
std::vector< Number * > inv_jacobian
unsigned int local_q_point_id_host(const unsigned int cell, const DataHost< dim, Number > &data, const unsigned int n_q_points, const unsigned int i)
void cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
Abstract base class for mapping classes.
Definition: mapping.h:311
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:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
constexpr int warp_size
Definition: cuda_size.h:40
constexpr __host__ unsigned int cells_per_block_shmem(int dim, int fe_degree)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void copy_to_host(const ArrayView< const T, MemorySpace::CUDA > &in, ArrayView< T, MemorySpace::Host > &out)
Definition: cuda.h:132
unsigned int global_dof_index
Definition: types.h:76
std::vector<::internal::MatrixFreeFunctions::ConstraintKinds > constraint_mask
std::vector< types::global_dof_index > local_to_global
std::vector< Point< dim, Number > > q_points
std::vector< Number > inv_jacobian
std::vector< Number > JxW
AdditionalData(const ParallelizationScheme parallelization_scheme=parallel_in_elem, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values|update_quadrature_points, const bool use_coloring=false, const bool overlap_communication_computation=false)
::internal::MatrixFreeFunctions::ConstraintKinds * constraint_mask
types::global_dof_index * local_to_global
SharedData(Number *vd, Number *gq[dim])
const MPI_Comm & comm