Reference documentation for deal.II version Git 44bda21415 2021-03-04 10:32:54 +0100
\(\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 - 2020 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 
43 
44 namespace CUDAWrappers
45 {
46  // forward declaration
47 # ifndef DOXYGEN
48  namespace internal
49  {
50  template <int dim, typename Number>
51  class ReinitHelper;
52  }
53 # endif
54 
81  template <int dim, typename Number = double>
82  class MatrixFree : public Subscriptor
83  {
84  public:
87  using CellFilter =
89 
96  {
98  parallel_over_elem
99  };
100 
105  {
110  const ParallelizationScheme parallelization_scheme = parallel_in_elem,
111  const UpdateFlags mapping_update_flags = update_gradients |
114  const bool use_coloring = false,
115  const bool overlap_communication_computation = false)
116  : parallelization_scheme(parallelization_scheme)
117  , mapping_update_flags(mapping_update_flags)
118  , use_coloring(use_coloring)
119  , overlap_communication_computation(overlap_communication_computation)
120  {
121 # ifndef DEAL_II_MPI_WITH_CUDA_SUPPORT
122  AssertThrow(
123  overlap_communication_computation == false,
124  ExcMessage(
125  "Overlapping communication and computation requires CUDA-aware MPI."));
126 # endif
127  if (overlap_communication_computation == true)
128  AssertThrow(
129  use_coloring == false || overlap_communication_computation == false,
130  ExcMessage(
131  "Overlapping communication and coloring are incompatible options. Only one of them can be enabled."));
132  }
148 
155 
161  };
162 
167  struct Data
168  {
173 
179 
183  Number *inv_jacobian;
184 
188  Number *JxW;
189 
193  unsigned int id;
194 
198  unsigned int n_cells;
199 
203  unsigned int padding_length;
204 
208  unsigned int row_start;
209 
213  unsigned int *constraint_mask;
214 
220  };
221 
225  MatrixFree();
226 
230  ~MatrixFree();
231 
235  unsigned int
236  get_padding_length() const;
237 
248  template <typename IteratorFiltersType>
249  void
250  reinit(const Mapping<dim> & mapping,
251  const DoFHandler<dim> & dof_handler,
252  const AffineConstraints<Number> &constraints,
253  const Quadrature<1> & quad,
254  const IteratorFiltersType & iterator_filter,
255  const AdditionalData & additional_data = AdditionalData());
256 
260  void
261  reinit(const Mapping<dim> & mapping,
262  const DoFHandler<dim> & dof_handler,
263  const AffineConstraints<Number> &constraints,
264  const Quadrature<1> & quad,
265  const AdditionalData & additional_data = AdditionalData());
266 
270  void
271  reinit(const DoFHandler<dim> & dof_handler,
272  const AffineConstraints<Number> &constraints,
273  const Quadrature<1> & quad,
275 
279  Data
280  get_data(unsigned int color) const;
281 
282  // clang-format off
300  // clang-format on
301  template <typename Functor, typename VectorType>
302  void
303  cell_loop(const Functor & func,
304  const VectorType &src,
305  VectorType & dst) const;
306 
322  template <typename Functor>
323  void
324  evaluate_coefficients(Functor func) const;
325 
330  template <typename VectorType>
331  void
332  copy_constrained_values(const VectorType &src, VectorType &dst) const;
333 
339  template <typename VectorType>
340  void
341  set_constrained_values(const Number val, VectorType &dst) const;
342 
347  void
348  initialize_dof_vector(
350 
356  void
357  initialize_dof_vector(
359 
363  const std::vector<std::vector<CellFilter>> &
364  get_colored_graph() const;
365 
376  const std::shared_ptr<const Utilities::MPI::Partitioner> &
377  get_vector_partitioner() const;
378 
382  void
383  free();
384 
388  const DoFHandler<dim> &
389  get_dof_handler() const;
390 
394  std::size_t
395  memory_consumption() const;
396 
397  private:
401  template <typename IteratorFiltersType>
402  void
403  internal_reinit(const Mapping<dim> & mapping,
404  const DoFHandler<dim> & dof_handler,
405  const AffineConstraints<Number> &constraints,
406  const Quadrature<1> & quad,
407  const IteratorFiltersType & iterator_filter,
408  std::shared_ptr<const MPI_Comm> comm,
409  const AdditionalData additional_data);
410 
415  template <typename Functor, typename VectorType>
416  void
417  serial_cell_loop(const Functor & func,
418  const VectorType &src,
419  VectorType & dst) const;
420 
425  template <typename Functor>
426  void
427  distributed_cell_loop(
428  const Functor & func,
431 
437  template <typename Functor>
438  void
439  distributed_cell_loop(
440  const Functor & func,
443 
448  template <typename VectorType>
449  void
450  serial_copy_constrained_values(const VectorType &src,
451  VectorType & dst) const;
452 
457  void
458  distributed_copy_constrained_values(
461 
468  void
469  distributed_copy_constrained_values(
472 
477  template <typename VectorType>
478  void
479  serial_set_constrained_values(const Number val, VectorType &dst) const;
480 
485  void
486  distributed_set_constrained_values(
487  const Number val,
489 
496  void
497  distributed_set_constrained_values(
498  const Number val,
500 
504  int my_id;
505 
511 
518 
524 
529 
533  unsigned int fe_degree;
534 
538  unsigned int dofs_per_cell;
539 
543  unsigned int n_constrained_dofs;
544 
548  unsigned int q_points_per_cell;
549 
553  unsigned int n_colors;
554 
558  std::vector<unsigned int> n_cells;
559 
564  std::vector<point_type *> q_points;
565 
570  std::vector<types::global_dof_index *> local_to_global;
571 
576  std::vector<Number *> inv_jacobian;
577 
582  std::vector<Number *> JxW;
583 
588 
592  std::vector<unsigned int *> constraint_mask;
593 
598  std::vector<dim3> grid_dim;
599 
604  std::vector<dim3> block_dim;
605 
610  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
611 
615  unsigned int cells_per_block;
616 
622 
628 
633  unsigned int padding_length;
634 
638  std::vector<unsigned int> row_start;
639 
644 
648  std::vector<std::vector<CellFilter>> graph;
649 
650  friend class internal::ReinitHelper<dim, Number>;
651  };
652 
653 
654 
655  // TODO find a better place to put these things
656 
660  template <int dim, typename Number>
661  struct SharedData
662  {
666  __device__
667  SharedData(Number *vd, Number *gq[dim])
668  : values(vd)
669  {
670  for (int d = 0; d < dim; ++d)
671  gradients[d] = gq[d];
672  }
673 
677  Number *values;
678 
684  Number *gradients[dim];
685  };
686 
687 
688 
689  // This function determines the number of cells per block, possibly at compile
690  // time (by virtue of being 'constexpr')
691  // TODO this function should be rewritten using meta-programming
692  __host__ __device__ constexpr unsigned int
693  cells_per_block_shmem(int dim, int fe_degree)
694  {
695  /* clang-format off */
696  // We are limiting the number of threads according to the
697  // following formulas:
698  // - in 2D: `threads = cells * (k+1)^d <= 4*CUDAWrappers::warp_size`
699  // - in 3D: `threads = cells * (k+1)^d <= 2*CUDAWrappers::warp_size`
700  return dim==2 ? (fe_degree==1 ? CUDAWrappers::warp_size : // 128
701  fe_degree==2 ? CUDAWrappers::warp_size/4 : // 72
702  fe_degree==3 ? CUDAWrappers::warp_size/8 : // 64
703  fe_degree==4 ? CUDAWrappers::warp_size/8 : // 100
704  1) :
705  dim==3 ? (fe_degree==1 ? CUDAWrappers::warp_size/4 : // 64
706  fe_degree==2 ? CUDAWrappers::warp_size/16 : // 54
707  1) : 1;
708  /* clang-format on */
709  }
710 
711 
712  /*----------------------- Helper functions ---------------------------------*/
718  template <int dim>
719  __device__ inline unsigned int
720  q_point_id_in_cell(const unsigned int n_q_points_1d)
721  {
722  return (dim == 1 ?
723  threadIdx.x % n_q_points_1d :
724  dim == 2 ?
725  threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
726  threadIdx.x % n_q_points_1d +
727  n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
728  }
729 
730 
731 
738  template <int dim, typename Number>
739  __device__ inline unsigned int
741  const unsigned int cell,
742  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
743  const unsigned int n_q_points_1d,
744  const unsigned int n_q_points)
745  {
746  return (data->row_start / data->padding_length + cell) * n_q_points +
747  q_point_id_in_cell<dim>(n_q_points_1d);
748  }
749 
750 
751 
757  template <int dim, typename Number>
758  __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
760  const unsigned int cell,
761  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
762  const unsigned int n_q_points_1d)
763  {
764  return *(data->q_points + data->padding_length * cell +
765  q_point_id_in_cell<dim>(n_q_points_1d));
766  }
767 
772  template <int dim, typename Number>
773  struct DataHost
774  {
778  std::vector<Point<dim, Number>> q_points;
779 
784  std::vector<types::global_dof_index> local_to_global;
785 
789  std::vector<Number> inv_jacobian;
790 
794  std::vector<Number> JxW;
795 
799  unsigned int id;
800 
804  unsigned int n_cells;
805 
809  unsigned int padding_length;
810 
814  unsigned int row_start;
815 
819  std::vector<unsigned int> constraint_mask;
820 
826  };
827 
828 
829 
836  template <int dim, typename Number>
839  const typename ::CUDAWrappers::MatrixFree<dim, Number>::Data &data,
840  const UpdateFlags &update_flags)
841  {
842  DataHost<dim, Number> data_host;
843 
844  data_host.id = data.id;
845  data_host.n_cells = data.n_cells;
846  data_host.padding_length = data.padding_length;
847  data_host.row_start = data.row_start;
848  data_host.use_coloring = data.use_coloring;
849 
850  const unsigned int n_elements =
851  data_host.n_cells * data_host.padding_length;
852  if (update_flags & update_quadrature_points)
853  {
854  data_host.q_points.resize(n_elements);
855  Utilities::CUDA::copy_to_host(data.q_points, data_host.q_points);
856  }
857 
858  data_host.local_to_global.resize(n_elements);
859  Utilities::CUDA::copy_to_host(data.local_to_global,
860  data_host.local_to_global);
861 
862  if (update_flags & update_gradients)
863  {
864  data_host.inv_jacobian.resize(n_elements * dim * dim);
865  Utilities::CUDA::copy_to_host(data.inv_jacobian,
866  data_host.inv_jacobian);
867  }
868 
869  if (update_flags & update_JxW_values)
870  {
871  data_host.JxW.resize(n_elements);
872  Utilities::CUDA::copy_to_host(data.JxW, data_host.JxW);
873  }
874 
875  data_host.constraint_mask.resize(data_host.n_cells);
876  Utilities::CUDA::copy_to_host(data.constraint_mask,
877  data_host.constraint_mask);
878 
879  return data_host;
880  }
881 
882 
883 
889  template <int dim, typename Number>
890  inline unsigned int
891  local_q_point_id_host(const unsigned int cell,
892  const DataHost<dim, Number> &data,
893  const unsigned int n_q_points,
894  const unsigned int i)
895  {
896  return (data.row_start / data.padding_length + cell) * n_q_points + i;
897  }
898 
899 
900 
908  template <int dim, typename Number>
909  inline Point<dim, Number>
910  get_quadrature_point_host(const unsigned int cell,
911  const DataHost<dim, Number> &data,
912  const unsigned int i)
913  {
914  return data.q_points[data.padding_length * cell + i];
915  }
916 
917 
918  /*----------------------- Inline functions ---------------------------------*/
919 
920 # ifndef DOXYGEN
921 
922  template <int dim, typename Number>
923  inline const std::vector<std::vector<
926  {
927  return graph;
928  }
929 
930 
931 
932  template <int dim, typename Number>
933  inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
935  {
936  return partitioner;
937  }
938 
939 
940 
941  template <int dim, typename Number>
942  inline const DoFHandler<dim> &
944  {
945  Assert(dof_handler != nullptr, ExcNotInitialized());
946 
947  return *dof_handler;
948  }
949 
950 # endif
951 
952 } // namespace CUDAWrappers
953 
955 
956 #endif
957 #endif
Transformed quadrature weights.
constexpr int warp_size
Definition: cuda_size.h:40
ParallelizationScheme parallelization_scheme
std::shared_ptr< const Utilities::MPI::Partitioner > partitioner
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
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)
__host__ constexpr unsigned int cells_per_block_shmem(int dim, int fe_degree)
Point< dim, Number > get_quadrature_point_host(const unsigned int cell, const DataHost< dim, Number > &data, const unsigned int i)
Transformed quadrature points.
static ::ExceptionBase & ExcNotInitialized()
types::global_dof_index * constrained_dofs
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
void copy_to_host(const ArrayView< const T, MemorySpace::CUDA > &in, ArrayView< T, MemorySpace::Host > &out)
Definition: cuda.h:132
types::global_dof_index n_dofs
Definition: point.h:110
const std::vector< std::vector< CellFilter > > & get_colored_graph() const
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)
std::vector< std::vector< CellFilter > > graph
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)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< point_type * > q_points
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner() const
unsigned int q_point_id_in_cell(const unsigned int n_q_points_1d)
SharedData(Number *vd, Number *gq[dim])
#define Assert(cond, exc)
Definition: exceptions.h:1466
std::vector< unsigned int > row_start
UpdateFlags
std::vector< Number * > inv_jacobian
Point< dim, Number > point_type
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:394
std::vector< Number > inv_jacobian
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)
std::vector< dim3 > block_dim
std::vector< unsigned int * > constraint_mask
std::vector< types::global_dof_index > local_to_global
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< unsigned int > constraint_mask
std::vector< dim3 > grid_dim
std::vector< types::global_dof_index * > local_to_global
*braid_SplitCommworld & comm
Definition: tensor.h:449
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:393
Shape function gradients.
const DoFHandler< dim > & get_dof_handler() const
std::vector< Number * > JxW
void free(T *&pointer)
Definition: cuda.h:97
std::vector< unsigned int > n_cells
DataHost< dim, Number > copy_mf_data_to_host(const typename::CUDAWrappers::MatrixFree< dim, Number >::Data &data, const UpdateFlags &update_flags)
std::vector< Number > JxW
types::global_dof_index * local_to_global
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
const DoFHandler< dim > * dof_handler
std::vector< Point< dim, Number > > q_points