Reference documentation for deal.II version Git f15f581df6 2020-07-10 15:30:09 -0400
\(\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 # include <deal.II/lac/cuda_vector.h>
38 
39 
41 
42 namespace CUDAWrappers
43 {
44  // forward declaration
45 # ifndef DOXYGEN
46  namespace internal
47  {
48  template <int dim, typename Number>
49  class ReinitHelper;
50  }
51 # endif
52 
79  template <int dim, typename Number = double>
80  class MatrixFree : public Subscriptor
81  {
82  public:
85 
92  {
94  parallel_over_elem
95  };
96 
101  {
106  const ParallelizationScheme parallelization_scheme = parallel_in_elem,
107  const UpdateFlags mapping_update_flags = update_gradients |
110  const bool use_coloring = false,
111  const bool overlap_communication_computation = false)
112  : parallelization_scheme(parallelization_scheme)
113  , mapping_update_flags(mapping_update_flags)
114  , use_coloring(use_coloring)
115  , overlap_communication_computation(overlap_communication_computation)
116  {
117 # ifndef DEAL_II_MPI_WITH_CUDA_SUPPORT
118  AssertThrow(
119  overlap_communication_computation == false,
120  ExcMessage(
121  "Overlapping communication and computation requires CUDA-aware MPI."));
122 # endif
123  if (overlap_communication_computation == true)
124  AssertThrow(
125  use_coloring == false || overlap_communication_computation == false,
126  ExcMessage(
127  "Overlapping communication and coloring are incompatible options. Only one of them can be enabled."));
128  }
144 
151 
157  };
158 
163  struct Data
164  {
169 
175 
179  Number *inv_jacobian;
180 
184  Number *JxW;
185 
189  unsigned int id;
190 
194  unsigned int n_cells;
195 
199  unsigned int padding_length;
200 
204  unsigned int row_start;
205 
209  unsigned int *constraint_mask;
210 
216  };
217 
221  MatrixFree();
222 
226  ~MatrixFree();
227 
231  unsigned int
232  get_padding_length() const;
233 
242  void
243  reinit(const Mapping<dim> & mapping,
244  const DoFHandler<dim> & dof_handler,
245  const AffineConstraints<Number> &constraints,
246  const Quadrature<1> & quad,
247  const AdditionalData & additional_data = AdditionalData());
248 
252  void
253  reinit(const DoFHandler<dim> & dof_handler,
254  const AffineConstraints<Number> &constraints,
255  const Quadrature<1> & quad,
257 
261  Data
262  get_data(unsigned int color) const;
263 
264  // clang-format off
282  // clang-format on
283  template <typename Functor, typename VectorType>
284  void
285  cell_loop(const Functor & func,
286  const VectorType &src,
287  VectorType & dst) const;
288 
304  template <typename Functor>
305  void
306  evaluate_coefficients(Functor func) const;
307 
312  template <typename VectorType>
313  void
314  copy_constrained_values(const VectorType &src, VectorType &dst) const;
315 
321  template <typename VectorType>
322  void
323  set_constrained_values(const Number val, VectorType &dst) const;
324 
329  void
330  initialize_dof_vector(
332 
338  void
339  initialize_dof_vector(
341 
352  const std::shared_ptr<const Utilities::MPI::Partitioner> &
353  get_vector_partitioner() const;
354 
358  void
359  free();
360 
364  const DoFHandler<dim> &
365  get_dof_handler() const;
366 
370  std::size_t
371  memory_consumption() const;
372 
373  private:
377  void
378  internal_reinit(const Mapping<dim> & mapping,
379  const DoFHandler<dim> & dof_handler,
380  const AffineConstraints<Number> &constraints,
381  const Quadrature<1> & quad,
382  std::shared_ptr<const MPI_Comm> comm,
383  const AdditionalData additional_data);
384 
389  template <typename Functor, typename VectorType>
390  void
391  serial_cell_loop(const Functor & func,
392  const VectorType &src,
393  VectorType & dst) const;
394 
399  template <typename Functor>
400  void
401  distributed_cell_loop(
402  const Functor & func,
405 
411  template <typename Functor>
412  void
413  distributed_cell_loop(
414  const Functor & func,
417 
422  template <typename VectorType>
423  void
424  serial_copy_constrained_values(const VectorType &src,
425  VectorType & dst) const;
426 
431  void
432  distributed_copy_constrained_values(
435 
442  void
443  distributed_copy_constrained_values(
446 
451  template <typename VectorType>
452  void
453  serial_set_constrained_values(const Number val, VectorType &dst) const;
454 
459  void
460  distributed_set_constrained_values(
461  const Number val,
463 
470  void
471  distributed_set_constrained_values(
472  const Number val,
474 
478  int my_id;
479 
485 
492 
498 
503 
507  unsigned int fe_degree;
508 
512  unsigned int dofs_per_cell;
513 
517  unsigned int n_constrained_dofs;
518 
522  unsigned int q_points_per_cell;
523 
527  unsigned int n_colors;
528 
532  std::vector<unsigned int> n_cells;
533 
538  std::vector<point_type *> q_points;
539 
544  std::vector<types::global_dof_index *> local_to_global;
545 
550  std::vector<Number *> inv_jacobian;
551 
556  std::vector<Number *> JxW;
557 
562 
566  std::vector<unsigned int *> constraint_mask;
567 
572  std::vector<dim3> grid_dim;
573 
578  std::vector<dim3> block_dim;
579 
584  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
585 
589  unsigned int cells_per_block;
590 
596 
602 
607  unsigned int padding_length;
608 
612  std::vector<unsigned int> row_start;
613 
618 
619  friend class internal::ReinitHelper<dim, Number>;
620  };
621 
622 
623 
624  // TODO find a better place to put these things
625 
629  template <int dim, typename Number>
630  struct SharedData
631  {
635  __device__
636  SharedData(Number *vd, Number *gq[dim])
637  : values(vd)
638  {
639  for (int d = 0; d < dim; ++d)
640  gradients[d] = gq[d];
641  }
642 
646  Number *values;
647 
653  Number *gradients[dim];
654  };
655 
656 
657 
658  // This function determines the number of cells per block, possibly at compile
659  // time (by virtue of being 'constexpr')
660  // TODO this function should be rewritten using meta-programming
661  __host__ __device__ constexpr unsigned int
662  cells_per_block_shmem(int dim, int fe_degree)
663  {
664  /* clang-format off */
665  // We are limiting the number of threads according to the
666  // following formulas:
667  // - in 2D: `threads = cells * (k+1)^d <= 4*CUDAWrappers::warp_size`
668  // - in 3D: `threads = cells * (k+1)^d <= 2*CUDAWrappers::warp_size`
669  return dim==2 ? (fe_degree==1 ? CUDAWrappers::warp_size : // 128
670  fe_degree==2 ? CUDAWrappers::warp_size/4 : // 72
671  fe_degree==3 ? CUDAWrappers::warp_size/8 : // 64
672  fe_degree==4 ? CUDAWrappers::warp_size/8 : // 100
673  1) :
674  dim==3 ? (fe_degree==1 ? CUDAWrappers::warp_size/4 : // 64
675  fe_degree==2 ? CUDAWrappers::warp_size/16 : // 54
676  1) : 1;
677  /* clang-format on */
678  }
679 
680 
681  /*----------------------- Helper functions ---------------------------------*/
687  template <int dim>
688  __device__ inline unsigned int
689  q_point_id_in_cell(const unsigned int n_q_points_1d)
690  {
691  return (dim == 1 ?
692  threadIdx.x % n_q_points_1d :
693  dim == 2 ?
694  threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
695  threadIdx.x % n_q_points_1d +
696  n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
697  }
698 
699 
700 
707  template <int dim, typename Number>
708  __device__ inline unsigned int
710  const unsigned int cell,
711  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
712  const unsigned int n_q_points_1d,
713  const unsigned int n_q_points)
714  {
715  return (data->row_start / data->padding_length + cell) * n_q_points +
716  q_point_id_in_cell<dim>(n_q_points_1d);
717  }
718 
719 
720 
726  template <int dim, typename Number>
727  __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
729  const unsigned int cell,
730  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
731  const unsigned int n_q_points_1d)
732  {
733  return *(data->q_points + data->padding_length * cell +
734  q_point_id_in_cell<dim>(n_q_points_1d));
735  }
736 
737 
738  /*----------------------- Inline functions ---------------------------------*/
739 
740 # ifndef DOXYGEN
741 
742  template <int dim, typename Number>
743  const std::shared_ptr<const Utilities::MPI::Partitioner> &
745  {
746  return partitioner;
747  }
748 
749  template <int dim, typename Number>
750  const DoFHandler<dim> &
752  {
753  Assert(dof_handler != nullptr, ExcNotInitialized());
754 
755  return *dof_handler;
756  }
757 
758 # endif
759 
760 } // namespace CUDAWrappers
761 
763 
764 #endif
765 
766 #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)
Transformed quadrature points.
static ::ExceptionBase & ExcNotInitialized()
types::global_dof_index * constrained_dofs
#define AssertThrow(cond, exc)
Definition: exceptions.h:1513
types::global_dof_index n_dofs
Definition: point.h:110
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)
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:1403
std::vector< unsigned int > row_start
UpdateFlags
std::vector< Number * > inv_jacobian
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< dim3 > grid_dim
std::vector< types::global_dof_index * > local_to_global
*braid_SplitCommworld & comm
Definition: tensor.h:448
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
Shape function gradients.
const DoFHandler< dim > & get_dof_handler() const
std::vector< Number * > JxW
void free(T *&pointer)
Definition: cuda.h:96
std::vector< unsigned int > n_cells
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