Reference documentation for deal.II version 9.4.0
\(\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 - 2022 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
25# include <deal.II/base/mpi.h>
27# include <deal.II/base/tensor.h>
28
30
32# include <deal.II/fe/mapping.h>
34
36
40
41
42
44
45// Forward declaration
46namespace internal
47{
48 namespace MatrixFreeFunctions
49 {
50 enum class ConstraintKinds : std::uint16_t;
51 }
52} // namespace internal
53
54namespace CUDAWrappers
55{
56 // forward declaration
57# ifndef DOXYGEN
58 namespace internal
59 {
60 template <int dim, typename Number>
61 class ReinitHelper;
62 }
63# endif
64
91 template <int dim, typename Number = double>
92 class MatrixFree : public Subscriptor
93 {
94 public:
97 using CellFilter =
99
106 {
109 };
110
115 {
124 const bool use_coloring = false,
125 const bool overlap_communication_computation = false)
130 {
131# ifndef DEAL_II_MPI_WITH_CUDA_SUPPORT
135 "Overlapping communication and computation requires CUDA-aware MPI."));
136# endif
141 "Overlapping communication and coloring are incompatible options. Only one of them can be enabled."));
142 }
158
165
171 };
172
177 struct Data
178 {
183
189
194
198 Number *JxW;
199
203 unsigned int id;
204
208 unsigned int n_cells;
209
213 unsigned int padding_length;
214
218 unsigned int row_start;
219
224
230 };
231
236
241
245 unsigned int
247
258 template <typename IteratorFiltersType>
259 void
260 reinit(const Mapping<dim> & mapping,
262 const AffineConstraints<Number> &constraints,
263 const Quadrature<1> & quad,
264 const IteratorFiltersType & iterator_filter,
265 const AdditionalData & additional_data = AdditionalData());
266
270 void
271 reinit(const Mapping<dim> & mapping,
273 const AffineConstraints<Number> &constraints,
274 const Quadrature<1> & quad,
275 const AdditionalData & additional_data = AdditionalData());
276
280 void
282 const AffineConstraints<Number> &constraints,
283 const Quadrature<1> & quad,
284 const AdditionalData & additional_data = AdditionalData());
285
289 Data
290 get_data(unsigned int color) const;
291
292 // clang-format off
310 // clang-format on
311 template <typename Functor, typename VectorType>
312 void
313 cell_loop(const Functor & func,
314 const VectorType &src,
315 VectorType & dst) const;
316
332 template <typename Functor>
333 void
334 evaluate_coefficients(Functor func) const;
335
340 template <typename VectorType>
341 void
342 copy_constrained_values(const VectorType &src, VectorType &dst) const;
343
349 template <typename VectorType>
350 void
351 set_constrained_values(const Number val, VectorType &dst) const;
352
357 void
360
366 void
369
373 const std::vector<std::vector<CellFilter>> &
375
386 const std::shared_ptr<const Utilities::MPI::Partitioner> &
388
392 void
394
398 const DoFHandler<dim> &
400
404 std::size_t
406
407 private:
411 template <typename IteratorFiltersType>
412 void
415 const AffineConstraints<Number> &constraints,
416 const Quadrature<1> & quad,
417 const IteratorFiltersType & iterator_filter,
418 std::shared_ptr<const MPI_Comm> comm,
419 const AdditionalData additional_data);
420
425 template <typename Functor, typename VectorType>
426 void
427 serial_cell_loop(const Functor & func,
428 const VectorType &src,
429 VectorType & dst) const;
430
435 template <typename Functor>
436 void
438 const Functor & func,
441
447 template <typename Functor>
448 void
450 const Functor & func,
453
458 template <typename VectorType>
459 void
460 serial_copy_constrained_values(const VectorType &src,
461 VectorType & dst) const;
462
467 void
471
478 void
482
487 template <typename VectorType>
488 void
489 serial_set_constrained_values(const Number val, VectorType &dst) const;
490
495 void
497 const Number val,
499
506 void
508 const Number val,
510
514 int my_id;
515
521
528
534
539
543 unsigned int fe_degree;
544
548 unsigned int dofs_per_cell;
549
553 unsigned int n_constrained_dofs;
554
558 unsigned int q_points_per_cell;
559
563 unsigned int n_colors;
564
568 std::vector<unsigned int> n_cells;
569
574 std::vector<point_type *> q_points;
575
580 std::vector<types::global_dof_index *> local_to_global;
581
586 std::vector<Number *> inv_jacobian;
587
592 std::vector<Number *> JxW;
593
598
602 std::vector<::internal::MatrixFreeFunctions::ConstraintKinds *>
604
609 std::vector<dim3> grid_dim;
610
615 std::vector<dim3> block_dim;
616
621 std::shared_ptr<const Utilities::MPI::Partitioner> partitioner;
622
626 unsigned int cells_per_block;
627
633
639
644 unsigned int padding_length;
645
649 std::vector<unsigned int> row_start;
650
655
659 std::vector<std::vector<CellFilter>> graph;
660
661 friend class internal::ReinitHelper<dim, Number>;
662 };
663
664
665
666 // TODO find a better place to put these things
667
671 template <int dim, typename Number>
673 {
677 __device__
678 SharedData(Number *vd, Number *gq[dim])
679 : values(vd)
680 {
681 for (unsigned int d = 0; d < dim; ++d)
682 gradients[d] = gq[d];
683 }
684
688 Number *values;
689
695 Number *gradients[dim];
696 };
697
698
699
700 // This function determines the number of cells per block, possibly at compile
701 // time (by virtue of being 'constexpr')
702 // TODO this function should be rewritten using meta-programming
703 __host__ __device__ constexpr unsigned int
704 cells_per_block_shmem(int dim, int fe_degree)
705 {
706 /* clang-format off */
707 // We are limiting the number of threads according to the
708 // following formulas:
709 // - in 2D: `threads = cells * (k+1)^d <= 4*CUDAWrappers::warp_size`
710 // - in 3D: `threads = cells * (k+1)^d <= 2*CUDAWrappers::warp_size`
711 return dim==2 ? (fe_degree==1 ? CUDAWrappers::warp_size : // 128
712 fe_degree==2 ? CUDAWrappers::warp_size/4 : // 72
713 fe_degree==3 ? CUDAWrappers::warp_size/8 : // 64
714 fe_degree==4 ? CUDAWrappers::warp_size/8 : // 100
715 1) :
716 dim==3 ? (fe_degree==1 ? CUDAWrappers::warp_size/4 : // 64
717 fe_degree==2 ? CUDAWrappers::warp_size/16 : // 54
718 1) : 1;
719 /* clang-format on */
720 }
721
722
723 /*----------------------- Helper functions ---------------------------------*/
729 template <int dim>
730 __device__ inline unsigned int
731 q_point_id_in_cell(const unsigned int n_q_points_1d)
732 {
733 return (
734 dim == 1 ? threadIdx.x % n_q_points_1d :
735 dim == 2 ? threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
736 threadIdx.x % n_q_points_1d +
737 n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
738 }
739
740
741
748 template <int dim, typename Number>
749 __device__ inline unsigned int
751 const unsigned int cell,
753 const unsigned int n_q_points_1d,
754 const unsigned int n_q_points)
755 {
756 return (data->row_start / data->padding_length + cell) * n_q_points +
757 q_point_id_in_cell<dim>(n_q_points_1d);
758 }
759
760
761
767 template <int dim, typename Number>
768 __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
770 const unsigned int cell,
772 const unsigned int n_q_points_1d)
773 {
774 return *(data->q_points + data->padding_length * cell +
775 q_point_id_in_cell<dim>(n_q_points_1d));
776 }
777
782 template <int dim, typename Number>
783 struct DataHost
784 {
788 std::vector<Point<dim, Number>> q_points;
789
794 std::vector<types::global_dof_index> local_to_global;
795
799 std::vector<Number> inv_jacobian;
800
804 std::vector<Number> JxW;
805
809 unsigned int id;
810
814 unsigned int n_cells;
815
819 unsigned int padding_length;
820
824 unsigned int row_start;
825
829 std::vector<::internal::MatrixFreeFunctions::ConstraintKinds>
831
837 };
838
839
840
847 template <int dim, typename Number>
850 const typename ::CUDAWrappers::MatrixFree<dim, Number>::Data &data,
851 const UpdateFlags &update_flags)
852 {
853 DataHost<dim, Number> data_host;
854
855 data_host.id = data.id;
856 data_host.n_cells = data.n_cells;
857 data_host.padding_length = data.padding_length;
858 data_host.row_start = data.row_start;
859 data_host.use_coloring = data.use_coloring;
860
861 const unsigned int n_elements =
862 data_host.n_cells * data_host.padding_length;
863 if (update_flags & update_quadrature_points)
864 {
865 data_host.q_points.resize(n_elements);
866 Utilities::CUDA::copy_to_host(data.q_points, data_host.q_points);
867 }
868
869 data_host.local_to_global.resize(n_elements);
870 Utilities::CUDA::copy_to_host(data.local_to_global,
871 data_host.local_to_global);
872
873 if (update_flags & update_gradients)
874 {
875 data_host.inv_jacobian.resize(n_elements * dim * dim);
876 Utilities::CUDA::copy_to_host(data.inv_jacobian,
877 data_host.inv_jacobian);
878 }
879
880 if (update_flags & update_JxW_values)
881 {
882 data_host.JxW.resize(n_elements);
883 Utilities::CUDA::copy_to_host(data.JxW, data_host.JxW);
884 }
885
886 data_host.constraint_mask.resize(data_host.n_cells);
887 Utilities::CUDA::copy_to_host(data.constraint_mask,
888 data_host.constraint_mask);
889
890 return data_host;
891 }
892
893
894
900 template <int dim, typename Number>
901 inline unsigned int
902 local_q_point_id_host(const unsigned int cell,
903 const DataHost<dim, Number> &data,
904 const unsigned int n_q_points,
905 const unsigned int i)
906 {
907 return (data.row_start / data.padding_length + cell) * n_q_points + i;
908 }
909
910
911
919 template <int dim, typename Number>
920 inline Point<dim, Number>
921 get_quadrature_point_host(const unsigned int cell,
922 const DataHost<dim, Number> &data,
923 const unsigned int i)
924 {
925 return data.q_points[data.padding_length * cell + i];
926 }
927
928
929 /*----------------------- Inline functions ---------------------------------*/
930
931# ifndef DOXYGEN
932
933 template <int dim, typename Number>
934 inline const std::vector<std::vector<
937 {
938 return graph;
939 }
940
941
942
943 template <int dim, typename Number>
944 inline const std::shared_ptr<const Utilities::MPI::Partitioner> &
946 {
947 return partitioner;
948 }
949
950
951
952 template <int dim, typename Number>
953 inline const DoFHandler<dim> &
955 {
956 Assert(dof_handler != nullptr, ExcNotInitialized());
957
958 return *dof_handler;
959 }
960
961# endif
962
963} // namespace CUDAWrappers
964
966
967#endif
968#endif
void evaluate_coefficients(Functor func) const
void serial_copy_constrained_values(const VectorType &src, VectorType &dst) const
std::vector< Number * > JxW
const DoFHandler< dim > * dof_handler
void initialize_dof_vector(LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &vec) const
Data get_data(unsigned int color) const
unsigned int get_padding_length() const
void initialize_dof_vector(LinearAlgebra::CUDAWrappers::Vector< Number > &vec) const
void serial_set_constrained_values(const Number val, VectorType &dst) const
std::vector< dim3 > block_dim
void reinit(const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const AdditionalData &additional_data=AdditionalData())
void distributed_copy_constrained_values(const LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &src, LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &dst) const
std::vector< unsigned int > n_cells
ParallelizationScheme parallelization_scheme
void internal_reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const IteratorFiltersType &iterator_filter, std::shared_ptr< const MPI_Comm > comm, const AdditionalData additional_data)
const DoFHandler< dim > & get_dof_handler() const
const std::shared_ptr< const Utilities::MPI::Partitioner > & get_vector_partitioner() const
void reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const IteratorFiltersType &iterator_filter, const AdditionalData &additional_data=AdditionalData())
void serial_cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
std::size_t memory_consumption() const
std::vector<::internal::MatrixFreeFunctions::ConstraintKinds * > constraint_mask
std::vector< point_type * > q_points
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)
void distributed_set_constrained_values(const Number val, LinearAlgebra::CUDAWrappers::Vector< Number > &dst) const
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)
void distributed_set_constrained_values(const Number val, LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &dst) const
Point< dim, Number > get_quadrature_point_host(const unsigned int cell, const DataHost< dim, Number > &data, const unsigned int i)
void distributed_cell_loop(const Functor &func, const LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &src, LinearAlgebra::distributed::Vector< Number, MemorySpace::CUDA > &dst) const
const std::vector< std::vector< CellFilter > > & get_colored_graph() const
std::vector< dim3 > grid_dim
std::vector< types::global_dof_index * > local_to_global
std::vector< std::vector< CellFilter > > graph
void copy_constrained_values(const VectorType &src, VectorType &dst) const
std::shared_ptr< const Utilities::MPI::Partitioner > partitioner
unsigned int q_point_id_in_cell(const unsigned int n_q_points_1d)
void distributed_copy_constrained_values(const LinearAlgebra::CUDAWrappers::Vector< Number > &src, LinearAlgebra::CUDAWrappers::Vector< Number > &dst) const
void reinit(const Mapping< dim > &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< Number > &constraints, const Quadrature< 1 > &quad, const AdditionalData &additional_data=AdditionalData())
types::global_dof_index n_dofs
void distributed_cell_loop(const Functor &func, const LinearAlgebra::CUDAWrappers::Vector< Number > &src, LinearAlgebra::CUDAWrappers::Vector< Number > &dst) const
std::vector< unsigned int > row_start
void set_constrained_values(const Number val, VectorType &dst) const
DataHost< dim, Number > copy_mf_data_to_host(const typename::CUDAWrappers::MatrixFree< dim, Number >::Data &data, const UpdateFlags &update_flags)
types::global_dof_index * constrained_dofs
std::vector< Number * > inv_jacobian
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)
void cell_loop(const Functor &func, const VectorType &src, VectorType &dst) const
Abstract base class for mapping classes.
Definition: mapping.h:311
Definition: tensor.h:503
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
UpdateFlags
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
__host__ constexpr unsigned int cells_per_block_shmem(int dim, int fe_degree)
constexpr int warp_size
Definition: cuda_size.h:40
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void copy_to_host(const ArrayView< const T, MemorySpace::CUDA > &in, ArrayView< T, MemorySpace::Host > &out)
Definition: cuda.h:132
std::vector<::internal::MatrixFreeFunctions::ConstraintKinds > constraint_mask
std::vector< types::global_dof_index > local_to_global
std::vector< Point< dim, Number > > q_points
std::vector< Number > inv_jacobian
std::vector< Number > JxW
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)
::internal::MatrixFreeFunctions::ConstraintKinds * constraint_mask
types::global_dof_index * local_to_global
SharedData(Number *vd, Number *gq[dim])
const MPI_Comm & comm