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