Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
cuda_matrix_free.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2019 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/mpi.h>
25 # include <deal.II/base/quadrature.h>
26 # include <deal.II/base/tensor.h>
27 
28 # include <deal.II/dofs/dof_handler.h>
29 
30 # include <deal.II/fe/fe_update_flags.h>
31 # include <deal.II/fe/mapping.h>
32 # include <deal.II/fe/mapping_q1.h>
33 
34 # include <deal.II/lac/affine_constraints.h>
35 # include <deal.II/lac/cuda_vector.h>
36 # include <deal.II/lac/la_parallel_vector.h>
37 
38 
39 DEAL_II_NAMESPACE_OPEN
40 
41 namespace CUDAWrappers
42 {
43  // forward declaration
44  namespace internal
45  {
46  template <int dim, typename Number>
47  class ReinitHelper;
48  }
49 
76  template <int dim, typename Number = double>
77  class MatrixFree : public Subscriptor
78  {
79  public:
82 
89  {
90  parallel_in_elem,
91  parallel_over_elem
92  };
93 
94  struct AdditionalData
95  {
96  AdditionalData(
97  const ParallelizationScheme parallelization_scheme = parallel_in_elem,
98  const UpdateFlags mapping_update_flags = update_gradients |
100  : parallelization_scheme(parallelization_scheme)
101  , mapping_update_flags(mapping_update_flags)
102  {}
103 
107  unsigned int n_colors;
112  ParallelizationScheme parallelization_scheme;
122  UpdateFlags mapping_update_flags;
123  };
124 
129  struct Data
130  {
131  point_type * q_points;
132  types::global_dof_index *local_to_global;
133  Number * inv_jacobian;
134  Number * JxW;
135  unsigned int n_cells;
136  unsigned int padding_length;
137  unsigned int row_start;
138  unsigned int * constraint_mask;
139  };
140 
144  MatrixFree();
145 
149  unsigned int
150  get_padding_length() const;
151 
160  void
161  reinit(const Mapping<dim> & mapping,
162  const DoFHandler<dim> & dof_handler,
163  const AffineConstraints<Number> &constraints,
164  const Quadrature<1> & quad,
165  const AdditionalData additional_data = AdditionalData());
166 
170  void
171  reinit(const DoFHandler<dim> & dof_handler,
172  const AffineConstraints<Number> &constraints,
173  const Quadrature<1> & quad,
174  const AdditionalData AdditionalData = AdditionalData());
175 
179  Data
180  get_data(unsigned int color) const;
181 
182  // clang-format off
200  // clang-format on
201  template <typename Functor, typename VectorType>
202  void
203  cell_loop(const Functor & func,
204  const VectorType &src,
205  VectorType & dst) const;
206 
222  template <typename Functor>
223  void
224  evaluate_coefficients(Functor func) const;
225 
230  template <typename VectorType>
231  void
232  copy_constrained_values(const VectorType &src, VectorType &dst) const;
233 
239  template <typename VectorType>
240  void
241  set_constrained_values(const Number val, VectorType &dst) const;
242 
246  void
247  free();
248 
252  std::size_t
253  memory_consumption() const;
254 
255  private:
259  void
260  internal_reinit(const Mapping<dim> & mapping,
261  const DoFHandler<dim> & dof_handler,
262  const AffineConstraints<Number> &constraints,
263  const Quadrature<1> & quad,
264  std::shared_ptr<const MPI_Comm> comm,
265  const AdditionalData additional_data);
266 
271  template <typename Functor, typename VectorType>
272  void
273  serial_cell_loop(const Functor & func,
274  const VectorType &src,
275  VectorType & dst) const;
276 
281  template <typename Functor>
282  void
283  distributed_cell_loop(
284  const Functor & func,
287 
293  template <typename Functor>
294  void
295  distributed_cell_loop(
296  const Functor & func,
299 
304  template <typename VectorType>
305  void
306  serial_copy_constrained_values(const VectorType &src,
307  VectorType & dst) const;
308 
313  void
314  distributed_copy_constrained_values(
317 
324  void
325  distributed_copy_constrained_values(
328 
333  template <typename VectorType>
334  void
335  serial_set_constrained_values(const Number val, VectorType &dst) const;
336 
341  void
342  distributed_set_constrained_values(
343  const Number val,
345 
352  void
353  distributed_set_constrained_values(
354  const Number val,
356 
362 
366  unsigned int fe_degree;
367 
371  unsigned int dofs_per_cell;
372 
376  unsigned int n_constrained_dofs;
377 
381  unsigned int q_points_per_cell;
382 
386  unsigned int n_colors;
387 
391  std::vector<unsigned int> n_cells;
392 
397  std::vector<point_type *> q_points;
398 
403  std::vector<types::global_dof_index *> local_to_global;
404 
409  std::vector<Number *> inv_jacobian;
410 
415  std::vector<Number *> JxW;
416 
417  // Pointer to the constrained degrees of freedom.
418  types::global_dof_index *constrained_dofs;
419 
420  // Mask deciding where constraints are set on a given cell.
421  std::vector<unsigned int *> constraint_mask;
422 
427  std::vector<dim3> grid_dim;
428 
433  std::vector<dim3> block_dim;
434 
439  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
440 
441  // Parallelization parameters
442  unsigned int cells_per_block;
443  dim3 constraint_grid_dim;
444  dim3 constraint_block_dim;
445 
446  unsigned int padding_length;
447  std::vector<unsigned int> row_start;
448 
449  friend class internal::ReinitHelper<dim, Number>;
450  };
451 
452 
453 
454  // TODO find a better place to put these things
455  // Structure to pass the shared memory into a general user function.
456  template <int dim, typename Number>
457  struct SharedData
458  {
459  __device__
460  SharedData(Number *vd, Number *gq[dim])
461  : values(vd)
462  {
463  for (int d = 0; d < dim; ++d)
464  gradients[d] = gq[d];
465  }
466 
467  Number *values;
468  Number *gradients[dim];
469  };
470 
471 
472 
473  // This function determines the number of cells per block, possibly at compile
474  // time (by virtue of being 'constexpr')
475  // TODO this function should be rewritten using meta-programming
476  __host__ __device__ constexpr unsigned int
477  cells_per_block_shmem(int dim, int fe_degree)
478  {
479  /* clang-format off */
480  return dim==2 ? (fe_degree==1 ? 32 :
481  fe_degree==2 ? 8 :
482  fe_degree==3 ? 4 :
483  fe_degree==4 ? 4 :
484  1) :
485  dim==3 ? (fe_degree==1 ? 8 :
486  fe_degree==2 ? 2 :
487  1) : 1;
488  /* clang-format on */
489  }
490 
491 
492 
498  template <int dim>
499  __device__ inline unsigned int
500  q_point_id_in_cell(const unsigned int n_q_points_1d)
501  {
502  return (dim == 1 ?
503  threadIdx.x % n_q_points_1d :
504  dim == 2 ?
505  threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
506  threadIdx.x % n_q_points_1d +
507  n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
508  }
509 
510 
511 
518  template <int dim, typename Number>
519  __device__ inline unsigned int
521  const unsigned int cell,
522  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
523  const unsigned int n_q_points_1d,
524  const unsigned int n_q_points)
525  {
526  return (data->row_start / data->padding_length + cell) * n_q_points +
527  q_point_id_in_cell<dim>(n_q_points_1d);
528  }
529 
530 
531 
537  template <int dim, typename Number>
538  __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
540  const unsigned int cell,
541  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
542  const unsigned int n_q_points_1d)
543  {
544  return *(data->q_points + data->padding_length * cell +
545  q_point_id_in_cell<dim>(n_q_points_1d));
546  }
547 } // namespace CUDAWrappers
548 
549 DEAL_II_NAMESPACE_CLOSE
550 
551 #endif
552 
553 #endif
Transformed quadrature weights.
ParallelizationScheme parallelization_scheme
std::shared_ptr< const Utilities::MPI::Partitioner > partitioner
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)
Definition: point.h:110
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)
std::vector< point_type * > q_points
UpdateFlags
std::vector< Number * > inv_jacobian
std::vector< dim3 > block_dim
unsigned int global_dof_index
Definition: types.h:89
std::vector< dim3 > grid_dim
std::vector< types::global_dof_index * > local_to_global
Definition: mpi.h:90
Shape function gradients.
std::vector< Number * > JxW
std::vector< unsigned int > n_cells
unsigned int q_point_id_in_cell(const unsigned int n_q_points_1d)