Reference documentation for deal.II version Git 840703ecf3 2019-10-14 17:15:20 -0400
\(\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 # ifndef DOXYGEN
45  namespace internal
46  {
47  template <int dim, typename Number>
48  class ReinitHelper;
49  }
50 # endif
51 
78  template <int dim, typename Number = double>
79  class MatrixFree : public Subscriptor
80  {
81  public:
84 
91  {
92  parallel_in_elem,
93  parallel_over_elem
94  };
95 
100  {
105  const ParallelizationScheme parallelization_scheme = parallel_in_elem,
106  const UpdateFlags mapping_update_flags = update_gradients |
108  const bool use_coloring = false,
109  const bool n_colors = 1)
110  : parallelization_scheme(parallelization_scheme)
111  , mapping_update_flags(mapping_update_flags)
112  , use_coloring(use_coloring)
113  , n_colors(n_colors)
114  {}
115 
119  unsigned int n_colors;
135 
142  };
143 
148  struct Data
149  {
154 
160 
164  Number *inv_jacobian;
165 
169  Number *JxW;
170 
174  unsigned int n_cells;
175 
179  unsigned int padding_length;
180 
184  unsigned int row_start;
185 
189  unsigned int *constraint_mask;
190 
196  };
197 
201  MatrixFree();
202 
206  unsigned int
207  get_padding_length() const;
208 
217  void
218  reinit(const Mapping<dim> & mapping,
219  const DoFHandler<dim> & dof_handler,
220  const AffineConstraints<Number> &constraints,
221  const Quadrature<1> & quad,
222  const AdditionalData & additional_data = AdditionalData());
223 
227  void
228  reinit(const DoFHandler<dim> & dof_handler,
229  const AffineConstraints<Number> &constraints,
230  const Quadrature<1> & quad,
232 
236  Data
237  get_data(unsigned int color) const;
238 
239  // clang-format off
257  // clang-format on
258  template <typename Functor, typename VectorType>
259  void
260  cell_loop(const Functor & func,
261  const VectorType &src,
262  VectorType & dst) const;
263 
279  template <typename Functor>
280  void
281  evaluate_coefficients(Functor func) const;
282 
287  template <typename VectorType>
288  void
289  copy_constrained_values(const VectorType &src, VectorType &dst) const;
290 
296  template <typename VectorType>
297  void
298  set_constrained_values(const Number val, VectorType &dst) const;
299 
304  void
305  initialize_dof_vector(
307 
313  void
314  initialize_dof_vector(
316 
320  void
321  free();
322 
326  std::size_t
327  memory_consumption() const;
328 
329  private:
333  void
334  internal_reinit(const Mapping<dim> & mapping,
335  const DoFHandler<dim> & dof_handler,
336  const AffineConstraints<Number> &constraints,
337  const Quadrature<1> & quad,
338  std::shared_ptr<const MPI_Comm> comm,
339  const AdditionalData additional_data);
340 
345  template <typename Functor, typename VectorType>
346  void
347  serial_cell_loop(const Functor & func,
348  const VectorType &src,
349  VectorType & dst) const;
350 
355  template <typename Functor>
356  void
357  distributed_cell_loop(
358  const Functor & func,
361 
367  template <typename Functor>
368  void
369  distributed_cell_loop(
370  const Functor & func,
373 
378  template <typename VectorType>
379  void
380  serial_copy_constrained_values(const VectorType &src,
381  VectorType & dst) const;
382 
387  void
388  distributed_copy_constrained_values(
391 
398  void
399  distributed_copy_constrained_values(
402 
407  template <typename VectorType>
408  void
409  serial_set_constrained_values(const Number val, VectorType &dst) const;
410 
415  void
416  distributed_set_constrained_values(
417  const Number val,
419 
426  void
427  distributed_set_constrained_values(
428  const Number val,
430 
436 
443 
448 
452  unsigned int fe_degree;
453 
457  unsigned int dofs_per_cell;
458 
462  unsigned int n_constrained_dofs;
463 
467  unsigned int q_points_per_cell;
468 
472  unsigned int n_colors;
473 
477  std::vector<unsigned int> n_cells;
478 
483  std::vector<point_type *> q_points;
484 
489  std::vector<types::global_dof_index *> local_to_global;
490 
495  std::vector<Number *> inv_jacobian;
496 
501  std::vector<Number *> JxW;
502 
507 
511  std::vector<unsigned int *> constraint_mask;
512 
517  std::vector<dim3> grid_dim;
518 
523  std::vector<dim3> block_dim;
524 
529  std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
530 
534  unsigned int cells_per_block;
535 
541 
547 
552  unsigned int padding_length;
553 
557  std::vector<unsigned int> row_start;
558 
559  friend class internal::ReinitHelper<dim, Number>;
560  };
561 
562 
563 
564  // TODO find a better place to put these things
565 
569  template <int dim, typename Number>
570  struct SharedData
571  {
575  __device__
576  SharedData(Number *vd, Number *gq[dim])
577  : values(vd)
578  {
579  for (int d = 0; d < dim; ++d)
580  gradients[d] = gq[d];
581  }
582 
586  Number *values;
587 
593  Number *gradients[dim];
594  };
595 
596 
597 
598  // This function determines the number of cells per block, possibly at compile
599  // time (by virtue of being 'constexpr')
600  // TODO this function should be rewritten using meta-programming
601  __host__ __device__ constexpr unsigned int
602  cells_per_block_shmem(int dim, int fe_degree)
603  {
604  /* clang-format off */
605  return dim==2 ? (fe_degree==1 ? 32 :
606  fe_degree==2 ? 8 :
607  fe_degree==3 ? 4 :
608  fe_degree==4 ? 4 :
609  1) :
610  dim==3 ? (fe_degree==1 ? 8 :
611  fe_degree==2 ? 2 :
612  1) : 1;
613  /* clang-format on */
614  }
615 
616 
617 
623  template <int dim>
624  __device__ inline unsigned int
625  q_point_id_in_cell(const unsigned int n_q_points_1d)
626  {
627  return (dim == 1 ?
628  threadIdx.x % n_q_points_1d :
629  dim == 2 ?
630  threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
631  threadIdx.x % n_q_points_1d +
632  n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
633  }
634 
635 
636 
643  template <int dim, typename Number>
644  __device__ inline unsigned int
646  const unsigned int cell,
647  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
648  const unsigned int n_q_points_1d,
649  const unsigned int n_q_points)
650  {
651  return (data->row_start / data->padding_length + cell) * n_q_points +
652  q_point_id_in_cell<dim>(n_q_points_1d);
653  }
654 
655 
656 
662  template <int dim, typename Number>
663  __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
665  const unsigned int cell,
666  const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
667  const unsigned int n_q_points_1d)
668  {
669  return *(data->q_points + data->padding_length * cell +
670  q_point_id_in_cell<dim>(n_q_points_1d));
671  }
672 } // namespace CUDAWrappers
673 
674 DEAL_II_NAMESPACE_CLOSE
675 
676 #endif
677 
678 #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)
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