Reference documentation for deal.II version Git a0b41b6d0f 2020-02-26 20:08:13 -0600
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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/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 
29 # include <deal.II/dofs/dof_handler.h>
30 
31 # include <deal.II/fe/fe_update_flags.h>
32 # include <deal.II/fe/mapping.h>
33 # include <deal.II/fe/mapping_q1.h>
34 
35 # include <deal.II/lac/affine_constraints.h>
36 # include <deal.II/lac/cuda_vector.h>
37 # include <deal.II/lac/la_parallel_vector.h>
38 
39 
40 DEAL_II_NAMESPACE_OPEN
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  {
93  parallel_in_elem,
94  parallel_over_elem
95  };
96 
101  {
106  const ParallelizationScheme parallelization_scheme = parallel_in_elem,
107  const UpdateFlags mapping_update_flags = update_gradients |
109  const bool use_coloring = false,
110  const bool n_colors = 1)
111  : parallelization_scheme(parallelization_scheme)
112  , mapping_update_flags(mapping_update_flags)
113  , use_coloring(use_coloring)
114  , n_colors(n_colors)
115  {}
116 
120  unsigned int n_colors;
136 
143  };
144 
149  struct Data
150  {
155 
161 
165  Number *inv_jacobian;
166 
170  Number *JxW;
171 
175  unsigned int n_cells;
176 
180  unsigned int padding_length;
181 
185  unsigned int row_start;
186 
190  unsigned int *constraint_mask;
191 
197  };
198 
202  MatrixFree();
203 
207  unsigned int
208  get_padding_length() const;
209 
218  void
219  reinit(const Mapping<dim> & mapping,
220  const DoFHandler<dim> & dof_handler,
221  const AffineConstraints<Number> &constraints,
222  const Quadrature<1> & quad,
223  const AdditionalData & additional_data = AdditionalData());
224 
228  void
229  reinit(const DoFHandler<dim> & dof_handler,
230  const AffineConstraints<Number> &constraints,
231  const Quadrature<1> & quad,
233 
237  Data
238  get_data(unsigned int color) const;
239 
240  // clang-format off
258  // clang-format on
259  template <typename Functor, typename VectorType>
260  void
261  cell_loop(const Functor & func,
262  const VectorType &src,
263  VectorType & dst) const;
264 
280  template <typename Functor>
281  void
282  evaluate_coefficients(Functor func) const;
283 
288  template <typename VectorType>
289  void
290  copy_constrained_values(const VectorType &src, VectorType &dst) const;
291 
297  template <typename VectorType>
298  void
299  set_constrained_values(const Number val, VectorType &dst) const;
300 
305  void
306  initialize_dof_vector(
308 
314  void
315  initialize_dof_vector(
317 
321  void
322  free();
323 
327  std::size_t
328  memory_consumption() const;
329 
330  private:
334  void
335  internal_reinit(const Mapping<dim> & mapping,
336  const DoFHandler<dim> & dof_handler,
337  const AffineConstraints<Number> &constraints,
338  const Quadrature<1> & quad,
339  std::shared_ptr<const MPI_Comm> comm,
340  const AdditionalData additional_data);
341 
346  template <typename Functor, typename VectorType>
347  void
348  serial_cell_loop(const Functor & func,
349  const VectorType &src,
350  VectorType & dst) const;
351 
356  template <typename Functor>
357  void
358  distributed_cell_loop(
359  const Functor & func,
362 
368  template <typename Functor>
369  void
370  distributed_cell_loop(
371  const Functor & func,
374 
379  template <typename VectorType>
380  void
381  serial_copy_constrained_values(const VectorType &src,
382  VectorType & dst) const;
383 
388  void
389  distributed_copy_constrained_values(
392 
399  void
400  distributed_copy_constrained_values(
403 
408  template <typename VectorType>
409  void
410  serial_set_constrained_values(const Number val, VectorType &dst) const;
411 
416  void
417  distributed_set_constrained_values(
418  const Number val,
420 
427  void
428  distributed_set_constrained_values(
429  const Number val,
431 
437 
444 
449 
453  unsigned int fe_degree;
454 
458  unsigned int dofs_per_cell;
459 
463  unsigned int n_constrained_dofs;
464 
468  unsigned int q_points_per_cell;
469 
473  unsigned int n_colors;
474 
478  std::vector<unsigned int> n_cells;
479 
484  std::vector<point_type *> q_points;
485 
490  std::vector<types::global_dof_index *> local_to_global;
491 
496  std::vector<Number *> inv_jacobian;
497 
502  std::vector<Number *> JxW;
503 
508 
512  std::vector<unsigned int *> constraint_mask;
513 
518  std::vector<dim3> grid_dim;
519 
524  std::vector<dim3> block_dim;
525 
530  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
531 
535  unsigned int cells_per_block;
536 
542 
548 
553  unsigned int padding_length;
554 
558  std::vector<unsigned int> row_start;
559 
560  friend class internal::ReinitHelper<dim, Number>;
561  };
562 
563 
564 
565  // TODO find a better place to put these things
566 
570  template <int dim, typename Number>
571  struct SharedData
572  {
576  __device__
577  SharedData(Number *vd, Number *gq[dim])
578  : values(vd)
579  {
580  for (int d = 0; d < dim; ++d)
581  gradients[d] = gq[d];
582  }
583 
587  Number *values;
588 
594  Number *gradients[dim];
595  };
596 
597 
598 
599  // This function determines the number of cells per block, possibly at compile
600  // time (by virtue of being 'constexpr')
601  // TODO this function should be rewritten using meta-programming
602  __host__ __device__ constexpr unsigned int
603  cells_per_block_shmem(int dim, int fe_degree)
604  {
605  /* clang-format off */
606  // We are limiting the number of threads according to the
607  // following formulas:
608  // - in 2D: `threads = cells * (k+1)^d <= 4*CUDAWrappers::warp_size`
609  // - in 3D: `threads = cells * (k+1)^d <= 2*CUDAWrappers::warp_size`
610  return dim==2 ? (fe_degree==1 ? CUDAWrappers::warp_size : // 128
611  fe_degree==2 ? CUDAWrappers::warp_size/4 : // 72
612  fe_degree==3 ? CUDAWrappers::warp_size/8 : // 64
613  fe_degree==4 ? CUDAWrappers::warp_size/8 : // 100
614  1) :
615  dim==3 ? (fe_degree==1 ? CUDAWrappers::warp_size/4 : // 64
616  fe_degree==2 ? CUDAWrappers::warp_size/16 : // 54
617  1) : 1;
618  /* clang-format on */
619  }
620 
621 
622 
628  template <int dim>
629  __device__ inline unsigned int
630  q_point_id_in_cell(const unsigned int n_q_points_1d)
631  {
632  return (dim == 1 ?
633  threadIdx.x % n_q_points_1d :
634  dim == 2 ?
635  threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
636  threadIdx.x % n_q_points_1d +
637  n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
638  }
639 
640 
641 
648  template <int dim, typename Number>
649  __device__ inline unsigned int
651  const unsigned int cell,
652  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
653  const unsigned int n_q_points_1d,
654  const unsigned int n_q_points)
655  {
656  return (data->row_start / data->padding_length + cell) * n_q_points +
657  q_point_id_in_cell<dim>(n_q_points_1d);
658  }
659 
660 
661 
667  template <int dim, typename Number>
668  __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
670  const unsigned int cell,
671  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
672  const unsigned int n_q_points_1d)
673  {
674  return *(data->q_points + data->padding_length * cell +
675  q_point_id_in_cell<dim>(n_q_points_1d));
676  }
677 } // namespace CUDAWrappers
678 
679 DEAL_II_NAMESPACE_CLOSE
680 
681 #endif
682 
683 #endif
Transformed quadrature weights.
constexpr int warp_size
Definition: cuda_size.h:40
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)
AdditionalData(const ParallelizationScheme parallelization_scheme=parallel_in_elem, const UpdateFlags mapping_update_flags=update_gradients|update_JxW_values, const bool use_coloring=false, const bool n_colors=1)
types::global_dof_index * constrained_dofs
types::global_dof_index n_dofs
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
SharedData(Number *vd, Number *gq[dim])
std::vector< unsigned int > row_start
UpdateFlags
std::vector< Number * > inv_jacobian
std::vector< dim3 > block_dim
std::vector< unsigned int * > constraint_mask
unsigned int global_dof_index
Definition: types.h:89
std::vector< dim3 > grid_dim
std::vector< types::global_dof_index * > local_to_global
Definition: tensor.h:422
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)
types::global_dof_index * local_to_global