Reference documentation for deal.II version 9.2.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_tensor_product_kernels.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 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_tensor_product_kernels_h
18 #define dealii_cuda_tensor_product_kernels_h
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/utilities.h>
23 
24 #include <deal.II/matrix_free/cuda_matrix_free.templates.h>
25 
26 #ifdef DEAL_II_COMPILER_CUDA_AWARE
27 
28 DEAL_II_NAMESPACE_OPEN
29 
30 
31 namespace CUDAWrappers
32 {
33  namespace internal
34  {
41  // TODO: for now only the general variant is implemented
42  enum EvaluatorVariant
43  {
44  evaluate_general,
45  evaluate_symmetric,
46  evaluate_evenodd
47  };
48 
49 
50 
56  template <EvaluatorVariant variant,
57  int dim,
58  int fe_degree,
59  int n_q_points_1d,
60  typename Number>
61  struct EvaluatorTensorProduct
62  {};
63 
64 
65 
72  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
73  struct EvaluatorTensorProduct<evaluate_general,
74  dim,
75  fe_degree,
76  n_q_points_1d,
77  Number>
78  {
79  static constexpr unsigned int dofs_per_cell =
80  Utilities::pow(fe_degree + 1, dim);
81  static constexpr unsigned int n_q_points =
82  Utilities::pow(n_q_points_1d, dim);
83 
84  __device__
85  EvaluatorTensorProduct();
86 
91  template <int direction, bool dof_to_quad, bool add, bool in_place>
92  __device__ void
93  values(Number shape_values[], const Number *in, Number *out) const;
94 
99  template <int direction, bool dof_to_quad, bool add, bool in_place>
100  __device__ void
101  gradients(Number shape_gradients[], const Number *in, Number *out) const;
102 
106  template <int direction, bool dof_to_quad, bool add, bool in_place>
107  __device__ void
108  apply(Number shape_data[], const Number *in, Number *out) const;
109 
113  __device__ void
114  value_at_quad_pts(Number *u);
115 
119  __device__ void
120  integrate_value(Number *u);
121 
126  __device__ void
127  gradient_at_quad_pts(const Number *const u, Number *grad_u[dim]);
128 
133  __device__ void
134  value_and_gradient_at_quad_pts(Number *const u, Number *grad_u[dim]);
135 
140  template <bool add>
141  __device__ void
142  integrate_gradient(Number *u, Number *grad_u[dim]);
143 
148  __device__ void
149  integrate_value_and_gradient(Number *u, Number *grad_u[dim]);
150  };
151 
152 
153 
154  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
155  __device__
156  EvaluatorTensorProduct<evaluate_general,
157  dim,
158  fe_degree,
159  n_q_points_1d,
160  Number>::EvaluatorTensorProduct()
161  {}
162 
163 
164 
165  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
166  template <int direction, bool dof_to_quad, bool add, bool in_place>
167  __device__ void
168  EvaluatorTensorProduct<evaluate_general,
169  dim,
170  fe_degree,
171  n_q_points_1d,
172  Number>::values(Number shape_values[],
173  const Number *in,
174  Number * out) const
175  {
176  apply<direction, dof_to_quad, add, in_place>(shape_values, in, out);
177  }
178 
179 
180 
181  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
182  template <int direction, bool dof_to_quad, bool add, bool in_place>
183  __device__ void
184  EvaluatorTensorProduct<evaluate_general,
185  dim,
186  fe_degree,
187  n_q_points_1d,
188  Number>::gradients(Number shape_gradients[],
189  const Number *in,
190  Number * out) const
191  {
192  apply<direction, dof_to_quad, add, in_place>(shape_gradients, in, out);
193  }
194 
195 
196 
197  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
198  template <int direction, bool dof_to_quad, bool add, bool in_place>
199  __device__ void
200  EvaluatorTensorProduct<evaluate_general,
201  dim,
202  fe_degree,
203  n_q_points_1d,
204  Number>::apply(Number shape_data[],
205  const Number *in,
206  Number * out) const
207  {
208  const unsigned int i = (dim == 1) ? 0 : threadIdx.x % n_q_points_1d;
209  const unsigned int j = (dim == 3) ? threadIdx.y : 0;
210  const unsigned int q = (dim == 1) ?
211  (threadIdx.x % n_q_points_1d) :
212  (dim == 2) ? threadIdx.y : threadIdx.z;
213 
214  // This loop simply multiply the shape function at the quadrature point by
215  // the value finite element coefficient.
216  Number t = 0;
217  for (int k = 0; k < n_q_points_1d; ++k)
218  {
219  const unsigned int shape_idx =
220  dof_to_quad ? (q + k * n_q_points_1d) : (k + q * n_q_points_1d);
221  const unsigned int source_idx =
222  (direction == 0) ?
223  (k + n_q_points_1d * (i + n_q_points_1d * j)) :
224  (direction == 1) ? (i + n_q_points_1d * (k + n_q_points_1d * j)) :
225  (i + n_q_points_1d * (j + n_q_points_1d * k));
226  t += shape_data[shape_idx] *
227  (in_place ? out[source_idx] : in[source_idx]);
228  }
229 
230  if (in_place)
231  __syncthreads();
232 
233  const unsigned int destination_idx =
234  (direction == 0) ?
235  (q + n_q_points_1d * (i + n_q_points_1d * j)) :
236  (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) :
237  (i + n_q_points_1d * (j + n_q_points_1d * q));
238 
239  if (add)
240  out[destination_idx] += t;
241  else
242  out[destination_idx] = t;
243  }
244 
245 
246 
247  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
248  inline __device__ void
249  EvaluatorTensorProduct<evaluate_general,
250  dim,
251  fe_degree,
252  n_q_points_1d,
253  Number>::value_at_quad_pts(Number *u)
254  {
255  switch (dim)
256  {
257  case 1:
258  {
259  values<0, true, false, true>(get_global_shape_values<Number>(),
260  u,
261  u);
262 
263  break;
264  }
265  case 2:
266  {
267  values<0, true, false, true>(get_global_shape_values<Number>(),
268  u,
269  u);
270  __syncthreads();
271  values<1, true, false, true>(get_global_shape_values<Number>(),
272  u,
273  u);
274 
275  break;
276  }
277  case 3:
278  {
279  values<0, true, false, true>(get_global_shape_values<Number>(),
280  u,
281  u);
282  __syncthreads();
283  values<1, true, false, true>(get_global_shape_values<Number>(),
284  u,
285  u);
286  __syncthreads();
287  values<2, true, false, true>(get_global_shape_values<Number>(),
288  u,
289  u);
290 
291  break;
292  }
293  default:
294  {
295  // Do nothing. We should throw but we can't from a __device__
296  // function.
297  }
298  }
299  }
300 
301 
302 
303  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
304  inline __device__ void
305  EvaluatorTensorProduct<evaluate_general,
306  dim,
307  fe_degree,
308  n_q_points_1d,
309  Number>::integrate_value(Number *u)
310  {
311  switch (dim)
312  {
313  case 1:
314  {
315  values<0, false, false, true>(get_global_shape_values<Number>(),
316  u,
317  u);
318 
319  break;
320  }
321  case 2:
322  {
323  values<0, false, false, true>(get_global_shape_values<Number>(),
324  u,
325  u);
326  __syncthreads();
327  values<1, false, false, true>(get_global_shape_values<Number>(),
328  u,
329  u);
330 
331  break;
332  }
333  case 3:
334  {
335  values<0, false, false, true>(get_global_shape_values<Number>(),
336  u,
337  u);
338  __syncthreads();
339  values<1, false, false, true>(get_global_shape_values<Number>(),
340  u,
341  u);
342  __syncthreads();
343  values<2, false, false, true>(get_global_shape_values<Number>(),
344  u,
345  u);
346 
347  break;
348  }
349  default:
350  {
351  // Do nothing. We should throw but we can't from a __device__
352  // function.
353  }
354  }
355  }
356 
357 
358 
359  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
360  inline __device__ void
361  EvaluatorTensorProduct<evaluate_general,
362  dim,
363  fe_degree,
364  n_q_points_1d,
365  Number>::gradient_at_quad_pts(const Number *const u,
366  Number *grad_u[dim])
367  {
368  switch (dim)
369  {
370  case 1:
371  {
372  gradients<0, true, false, false>(
373  get_global_shape_gradients<Number>(), u, grad_u[0]);
374 
375  break;
376  }
377  case 2:
378  {
379  gradients<0, true, false, false>(
380  get_global_shape_gradients<Number>(), u, grad_u[0]);
381  values<0, true, false, false>(get_global_shape_values<Number>(),
382  u,
383  grad_u[1]);
384 
385  __syncthreads();
386 
387  values<1, true, false, true>(get_global_shape_values<Number>(),
388  grad_u[0],
389  grad_u[0]);
390  gradients<1, true, false, true>(
391  get_global_shape_gradients<Number>(), grad_u[1], grad_u[1]);
392 
393  break;
394  }
395  case 3:
396  {
397  gradients<0, true, false, false>(
398  get_global_shape_gradients<Number>(), u, grad_u[0]);
399  values<0, true, false, false>(get_global_shape_values<Number>(),
400  u,
401  grad_u[1]);
402  values<0, true, false, false>(get_global_shape_values<Number>(),
403  u,
404  grad_u[2]);
405 
406  __syncthreads();
407 
408  values<1, true, false, true>(get_global_shape_values<Number>(),
409  grad_u[0],
410  grad_u[0]);
411  gradients<1, true, false, true>(
412  get_global_shape_gradients<Number>(), grad_u[1], grad_u[1]);
413  values<1, true, false, true>(get_global_shape_values<Number>(),
414  grad_u[2],
415  grad_u[2]);
416 
417  __syncthreads();
418 
419  values<2, true, false, true>(get_global_shape_values<Number>(),
420  grad_u[0],
421  grad_u[0]);
422  values<2, true, false, true>(get_global_shape_values<Number>(),
423  grad_u[1],
424  grad_u[1]);
425  gradients<2, true, false, true>(
426  get_global_shape_gradients<Number>(), grad_u[2], grad_u[2]);
427 
428  break;
429  }
430  default:
431  {
432  // Do nothing. We should throw but we can't from a __device__
433  // function.
434  }
435  }
436  }
437 
438 
439 
440  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
441  inline __device__ void
442  EvaluatorTensorProduct<
443  evaluate_general,
444  dim,
445  fe_degree,
446  n_q_points_1d,
447  Number>::value_and_gradient_at_quad_pts(Number *const u,
448  Number * grad_u[dim])
449  {
450  switch (dim)
451  {
452  case 1:
453  {
454  values<0, true, false, true>(get_global_shape_values<Number>(),
455  u,
456  u);
457  __syncthreads();
458 
459  gradients<0, true, false, false>(
460  get_global_co_shape_gradients<Number>(), u, grad_u[0]);
461 
462  break;
463  }
464  case 2:
465  {
466  values<0, true, false, true>(get_global_shape_values<Number>(),
467  u,
468  u);
469  __syncthreads();
470  values<1, true, false, true>(get_global_shape_values<Number>(),
471  u,
472  u);
473  __syncthreads();
474 
475  gradients<0, true, false, false>(
476  get_global_co_shape_gradients<Number>(), u, grad_u[0]);
477  gradients<1, true, false, false>(
478  get_global_co_shape_gradients<Number>(), u, grad_u[1]);
479 
480  break;
481  }
482  case 3:
483  {
484  values<0, true, false, true>(get_global_shape_values<Number>(),
485  u,
486  u);
487  __syncthreads();
488  values<1, true, false, true>(get_global_shape_values<Number>(),
489  u,
490  u);
491  __syncthreads();
492  values<2, true, false, true>(get_global_shape_values<Number>(),
493  u,
494  u);
495  __syncthreads();
496 
497  gradients<0, true, false, false>(
498  get_global_co_shape_gradients<Number>(), u, grad_u[0]);
499  gradients<1, true, false, false>(
500  get_global_co_shape_gradients<Number>(), u, grad_u[1]);
501  gradients<2, true, false, false>(
502  get_global_co_shape_gradients<Number>(), u, grad_u[2]);
503 
504  break;
505  }
506  default:
507  {
508  // Do nothing. We should throw but we can't from a __device__
509  // function.
510  }
511  }
512  }
513 
514 
515 
516  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
517  template <bool add>
518  inline __device__ void
519  EvaluatorTensorProduct<evaluate_general,
520  dim,
521  fe_degree,
522  n_q_points_1d,
523  Number>::integrate_gradient(Number *u,
524  Number *grad_u[dim])
525  {
526  switch (dim)
527  {
528  case 1:
529  {
530  gradients<0, false, add, false>(
531  get_global_shape_gradients<Number>(), grad_u[dim], u);
532 
533  break;
534  }
535  case 2:
536  {
537  gradients<0, false, false, true>(
538  get_global_shape_gradients<Number>(), grad_u[0], grad_u[0]);
539  values<0, false, false, true>(get_global_shape_values<Number>(),
540  grad_u[1],
541  grad_u[1]);
542 
543  __syncthreads();
544 
545  values<1, false, add, false>(get_global_shape_values<Number>(),
546  grad_u[0],
547  u);
548  __syncthreads();
549  gradients<1, false, true, false>(
550  get_global_shape_gradients<Number>(), grad_u[1], u);
551 
552  break;
553  }
554  case 3:
555  {
556  gradients<0, false, false, true>(
557  get_global_shape_gradients<Number>(), grad_u[0], grad_u[0]);
558  values<0, false, false, true>(get_global_shape_values<Number>(),
559  grad_u[1],
560  grad_u[1]);
561  values<0, false, false, true>(get_global_shape_values<Number>(),
562  grad_u[2],
563  grad_u[2]);
564 
565  __syncthreads();
566 
567  values<1, false, false, true>(get_global_shape_values<Number>(),
568  grad_u[0],
569  grad_u[0]);
570  gradients<1, false, false, true>(
571  get_global_shape_gradients<Number>(), grad_u[1], grad_u[1]);
572  values<1, false, false, true>(get_global_shape_values<Number>(),
573  grad_u[2],
574  grad_u[2]);
575 
576  __syncthreads();
577 
578  values<2, false, add, false>(get_global_shape_values<Number>(),
579  grad_u[0],
580  u);
581  __syncthreads();
582  values<2, false, true, false>(get_global_shape_values<Number>(),
583  grad_u[1],
584  u);
585  __syncthreads();
586  gradients<2, false, true, false>(
587  get_global_shape_gradients<Number>(), grad_u[2], u);
588 
589  break;
590  }
591  default:
592  {
593  // Do nothing. We should throw but we can't from a __device__
594  // function.
595  }
596  }
597  }
598 
599 
600 
601  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
602  inline __device__ void
603  EvaluatorTensorProduct<evaluate_general,
604  dim,
605  fe_degree,
606  n_q_points_1d,
607  Number>::integrate_value_and_gradient(Number *u,
608  Number
609  *grad_u[dim])
610  {
611  switch (dim)
612  {
613  case 1:
614  {
615  gradients<0, false, true, false>(
616  get_global_co_shape_gradients<Number>(), grad_u[0], u);
617  __syncthreads();
618 
619  values<0, false, false, true>(get_global_shape_values<Number>(),
620  u,
621  u);
622 
623  break;
624  }
625  case 2:
626  {
627  gradients<1, false, true, false>(
628  get_global_co_shape_gradients<Number>(), grad_u[1], u);
629  __syncthreads();
630  gradients<0, false, true, false>(
631  get_global_co_shape_gradients<Number>(), grad_u[0], u);
632  __syncthreads();
633 
634  values<1, false, false, true>(get_global_shape_values<Number>(),
635  u,
636  u);
637  __syncthreads();
638  values<0, false, false, true>(get_global_shape_values<Number>(),
639  u,
640  u);
641  __syncthreads();
642 
643  break;
644  }
645  case 3:
646  {
647  gradients<2, false, true, false>(
648  get_global_co_shape_gradients<Number>(), grad_u[2], u);
649  __syncthreads();
650  gradients<1, false, true, false>(
651  get_global_co_shape_gradients<Number>(), grad_u[1], u);
652  __syncthreads();
653  gradients<0, false, true, false>(
654  get_global_co_shape_gradients<Number>(), grad_u[0], u);
655  __syncthreads();
656 
657  values<2, false, false, true>(get_global_shape_values<Number>(),
658  u,
659  u);
660  __syncthreads();
661  values<1, false, false, true>(get_global_shape_values<Number>(),
662  u,
663  u);
664  __syncthreads();
665  values<0, false, false, true>(get_global_shape_values<Number>(),
666  u,
667  u);
668  __syncthreads();
669 
670  break;
671  }
672  default:
673  {
674  // Do nothing. We should throw but we can't from a __device__
675  // function.
676  }
677  }
678  }
679  } // namespace internal
680 } // namespace CUDAWrappers
681 
682 DEAL_II_NAMESPACE_CLOSE
683 
684 #endif
685 
686 #endif
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
CUDAWrappers::internal::evaluate_general
@ evaluate_general
Definition: cuda_tensor_product_kernels.h:44
utilities.h
CUDAWrappers::internal::EvaluatorTensorProduct
Definition: cuda_tensor_product_kernels.h:61
CUDAWrappers::internal::EvaluatorVariant
EvaluatorVariant
Definition: cuda_tensor_product_kernels.h:42
CUDAWrappers
Definition: cuda_size.h:23
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
CUDAWrappers::internal::evaluate_evenodd
@ evaluate_evenodd
Definition: cuda_tensor_product_kernels.h:46
CUDAWrappers::internal::evaluate_symmetric
@ evaluate_symmetric
Definition: cuda_tensor_product_kernels.h:45
config.h
internal
Definition: aligned_vector.h:369
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Generated by   doxygen 1.8.17