Reference documentation for deal.II version 9.3.3
\(\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 - 2020 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
28DEAL_II_NAMESPACE_OPEN
29
30
31namespace 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 const int mf_object_id;
64 };
65
66
67
74 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
75 struct EvaluatorTensorProduct<evaluate_general,
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__
87 EvaluatorTensorProduct(int mf_object_id);
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__
160 EvaluatorTensorProduct<evaluate_general,
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
173 EvaluatorTensorProduct<evaluate_general,
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
189 EvaluatorTensorProduct<evaluate_general,
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
205 EvaluatorTensorProduct<evaluate_general,
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) ?
216 (threadIdx.x % n_q_points_1d) :
217 (dim == 2) ? threadIdx.y : 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) ?
228 (k + n_q_points_1d * (i + n_q_points_1d * j)) :
229 (direction == 1) ? (i + n_q_points_1d * (k + n_q_points_1d * j)) :
230 (i + n_q_points_1d * (j + n_q_points_1d * k));
231 t += shape_data[shape_idx] *
232 (in_place ? out[source_idx] : in[source_idx]);
233 }
234
235 if (in_place)
236 __syncthreads();
237
238 const unsigned int destination_idx =
239 (direction == 0) ?
240 (q + n_q_points_1d * (i + n_q_points_1d * j)) :
241 (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) :
242 (i + n_q_points_1d * (j + n_q_points_1d * q));
243
244 if (add)
245 out[destination_idx] += t;
246 else
247 out[destination_idx] = t;
248 }
249
250
251
252 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
253 inline __device__ void
254 EvaluatorTensorProduct<evaluate_general,
255 dim,
256 fe_degree,
257 n_q_points_1d,
258 Number>::value_at_quad_pts(Number *u)
259 {
260 switch (dim)
261 {
262 case 1:
263 {
264 values<0, true, false, true>(
265 get_global_shape_values<Number>(mf_object_id), u, u);
266
267 break;
268 }
269 case 2:
270 {
271 values<0, true, false, true>(
272 get_global_shape_values<Number>(mf_object_id), u, u);
273 __syncthreads();
274 values<1, true, false, true>(
275 get_global_shape_values<Number>(mf_object_id), u, u);
276
277 break;
278 }
279 case 3:
280 {
281 values<0, true, false, true>(
282 get_global_shape_values<Number>(mf_object_id), u, u);
283 __syncthreads();
284 values<1, true, false, true>(
285 get_global_shape_values<Number>(mf_object_id), u, u);
286 __syncthreads();
287 values<2, true, false, true>(
288 get_global_shape_values<Number>(mf_object_id), u, u);
289
290 break;
291 }
292 default:
293 {
294 // Do nothing. We should throw but we can't from a __device__
295 // function.
296 }
297 }
298 }
299
300
301
302 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
303 inline __device__ void
304 EvaluatorTensorProduct<evaluate_general,
305 dim,
306 fe_degree,
307 n_q_points_1d,
308 Number>::integrate_value(Number *u)
309 {
310 switch (dim)
311 {
312 case 1:
313 {
314 values<0, false, false, true>(
315 get_global_shape_values<Number>(mf_object_id), u, u);
316
317 break;
318 }
319 case 2:
320 {
321 values<0, false, false, true>(
322 get_global_shape_values<Number>(mf_object_id), u, u);
323 __syncthreads();
324 values<1, false, false, true>(
325 get_global_shape_values<Number>(mf_object_id), u, u);
326
327 break;
328 }
329 case 3:
330 {
331 values<0, false, false, true>(
332 get_global_shape_values<Number>(mf_object_id), u, u);
333 __syncthreads();
334 values<1, false, false, true>(
335 get_global_shape_values<Number>(mf_object_id), u, u);
336 __syncthreads();
337 values<2, false, false, true>(
338 get_global_shape_values<Number>(mf_object_id), u, u);
339
340 break;
341 }
342 default:
343 {
344 // Do nothing. We should throw but we can't from a __device__
345 // function.
346 }
347 }
348 }
349
350
351
352 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
353 inline __device__ void
354 EvaluatorTensorProduct<evaluate_general,
355 dim,
356 fe_degree,
357 n_q_points_1d,
358 Number>::gradient_at_quad_pts(const Number *const u,
359 Number *grad_u[dim])
360 {
361 switch (dim)
362 {
363 case 1:
364 {
365 gradients<0, true, false, false>(
366 get_global_shape_gradients<Number>(mf_object_id), u, grad_u[0]);
367
368 break;
369 }
370 case 2:
371 {
372 gradients<0, true, false, false>(
373 get_global_shape_gradients<Number>(mf_object_id), u, grad_u[0]);
374 values<0, true, false, false>(
375 get_global_shape_values<Number>(mf_object_id), u, grad_u[1]);
376
377 __syncthreads();
378
379 values<1, true, false, true>(get_global_shape_values<Number>(
380 mf_object_id),
381 grad_u[0],
382 grad_u[0]);
383 gradients<1, true, false, true>(
384 get_global_shape_gradients<Number>(mf_object_id),
385 grad_u[1],
386 grad_u[1]);
387
388 break;
389 }
390 case 3:
391 {
392 gradients<0, true, false, false>(
393 get_global_shape_gradients<Number>(mf_object_id), u, grad_u[0]);
394 values<0, true, false, false>(
395 get_global_shape_values<Number>(mf_object_id), u, grad_u[1]);
396 values<0, true, false, false>(
397 get_global_shape_values<Number>(mf_object_id), u, grad_u[2]);
398
399 __syncthreads();
400
401 values<1, true, false, true>(get_global_shape_values<Number>(
402 mf_object_id),
403 grad_u[0],
404 grad_u[0]);
405 gradients<1, true, false, true>(
406 get_global_shape_gradients<Number>(mf_object_id),
407 grad_u[1],
408 grad_u[1]);
409 values<1, true, false, true>(get_global_shape_values<Number>(
410 mf_object_id),
411 grad_u[2],
412 grad_u[2]);
413
414 __syncthreads();
415
416 values<2, true, false, true>(get_global_shape_values<Number>(
417 mf_object_id),
418 grad_u[0],
419 grad_u[0]);
420 values<2, true, false, true>(get_global_shape_values<Number>(
421 mf_object_id),
422 grad_u[1],
423 grad_u[1]);
424 gradients<2, true, false, true>(
425 get_global_shape_gradients<Number>(mf_object_id),
426 grad_u[2],
427 grad_u[2]);
428
429 break;
430 }
431 default:
432 {
433 // Do nothing. We should throw but we can't from a __device__
434 // function.
435 }
436 }
437 }
438
439
440
441 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
442 inline __device__ void
443 EvaluatorTensorProduct<
444 evaluate_general,
445 dim,
446 fe_degree,
447 n_q_points_1d,
448 Number>::value_and_gradient_at_quad_pts(Number *const u,
449 Number * grad_u[dim])
450 {
451 switch (dim)
452 {
453 case 1:
454 {
455 values<0, true, false, true>(
456 get_global_shape_values<Number>(mf_object_id), u, u);
457 __syncthreads();
458
459 gradients<0, true, false, false>(
460 get_global_co_shape_gradients<Number>(mf_object_id),
461 u,
462 grad_u[0]);
463
464 break;
465 }
466 case 2:
467 {
468 values<0, true, false, true>(
469 get_global_shape_values<Number>(mf_object_id), u, u);
470 __syncthreads();
471 values<1, true, false, true>(
472 get_global_shape_values<Number>(mf_object_id), u, u);
473 __syncthreads();
474
475 gradients<0, true, false, false>(
476 get_global_co_shape_gradients<Number>(mf_object_id),
477 u,
478 grad_u[0]);
479 gradients<1, true, false, false>(
480 get_global_co_shape_gradients<Number>(mf_object_id),
481 u,
482 grad_u[1]);
483
484 break;
485 }
486 case 3:
487 {
488 values<0, true, false, true>(
489 get_global_shape_values<Number>(mf_object_id), u, u);
490 __syncthreads();
491 values<1, true, false, true>(
492 get_global_shape_values<Number>(mf_object_id), u, u);
493 __syncthreads();
494 values<2, true, false, true>(
495 get_global_shape_values<Number>(mf_object_id), u, u);
496 __syncthreads();
497
498 gradients<0, true, false, false>(
499 get_global_co_shape_gradients<Number>(mf_object_id),
500 u,
501 grad_u[0]);
502 gradients<1, true, false, false>(
503 get_global_co_shape_gradients<Number>(mf_object_id),
504 u,
505 grad_u[1]);
506 gradients<2, true, false, false>(
507 get_global_co_shape_gradients<Number>(mf_object_id),
508 u,
509 grad_u[2]);
510
511 break;
512 }
513 default:
514 {
515 // Do nothing. We should throw but we can't from a __device__
516 // function.
517 }
518 }
519 }
520
521
522
523 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
524 template <bool add>
525 inline __device__ void
526 EvaluatorTensorProduct<evaluate_general,
527 dim,
528 fe_degree,
529 n_q_points_1d,
530 Number>::integrate_gradient(Number *u,
531 Number *grad_u[dim])
532 {
533 switch (dim)
534 {
535 case 1:
536 {
537 gradients<0, false, add, false>(
538 get_global_shape_gradients<Number>(mf_object_id),
539 grad_u[dim],
540 u);
541
542 break;
543 }
544 case 2:
545 {
546 gradients<0, false, false, true>(
547 get_global_shape_gradients<Number>(mf_object_id),
548 grad_u[0],
549 grad_u[0]);
550 values<0, false, false, true>(get_global_shape_values<Number>(
551 mf_object_id),
552 grad_u[1],
553 grad_u[1]);
554
555 __syncthreads();
556
557 values<1, false, add, false>(
558 get_global_shape_values<Number>(mf_object_id), grad_u[0], u);
559 __syncthreads();
560 gradients<1, false, true, false>(
561 get_global_shape_gradients<Number>(mf_object_id), grad_u[1], u);
562
563 break;
564 }
565 case 3:
566 {
567 gradients<0, false, false, true>(
568 get_global_shape_gradients<Number>(mf_object_id),
569 grad_u[0],
570 grad_u[0]);
571 values<0, false, false, true>(get_global_shape_values<Number>(
572 mf_object_id),
573 grad_u[1],
574 grad_u[1]);
575 values<0, false, false, true>(get_global_shape_values<Number>(
576 mf_object_id),
577 grad_u[2],
578 grad_u[2]);
579
580 __syncthreads();
581
582 values<1, false, false, true>(get_global_shape_values<Number>(
583 mf_object_id),
584 grad_u[0],
585 grad_u[0]);
586 gradients<1, false, false, true>(
587 get_global_shape_gradients<Number>(mf_object_id),
588 grad_u[1],
589 grad_u[1]);
590 values<1, false, false, true>(get_global_shape_values<Number>(
591 mf_object_id),
592 grad_u[2],
593 grad_u[2]);
594
595 __syncthreads();
596
597 values<2, false, add, false>(
598 get_global_shape_values<Number>(mf_object_id), grad_u[0], u);
599 __syncthreads();
600 values<2, false, true, false>(
601 get_global_shape_values<Number>(mf_object_id), grad_u[1], u);
602 __syncthreads();
603 gradients<2, false, true, false>(
604 get_global_shape_gradients<Number>(mf_object_id), grad_u[2], u);
605
606 break;
607 }
608 default:
609 {
610 // Do nothing. We should throw but we can't from a __device__
611 // function.
612 }
613 }
614 }
615
616
617
618 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
619 inline __device__ void
620 EvaluatorTensorProduct<evaluate_general,
621 dim,
622 fe_degree,
623 n_q_points_1d,
624 Number>::integrate_value_and_gradient(Number *u,
625 Number
626 *grad_u[dim])
627 {
628 switch (dim)
629 {
630 case 1:
631 {
632 gradients<0, false, true, false>(
633 get_global_co_shape_gradients<Number>(mf_object_id),
634 grad_u[0],
635 u);
636 __syncthreads();
637
638 values<0, false, false, true>(
639 get_global_shape_values<Number>(mf_object_id), u, u);
640
641 break;
642 }
643 case 2:
644 {
645 gradients<1, false, true, false>(
646 get_global_co_shape_gradients<Number>(mf_object_id),
647 grad_u[1],
648 u);
649 __syncthreads();
650 gradients<0, false, true, false>(
651 get_global_co_shape_gradients<Number>(mf_object_id),
652 grad_u[0],
653 u);
654 __syncthreads();
655
656 values<1, false, false, true>(
657 get_global_shape_values<Number>(mf_object_id), u, u);
658 __syncthreads();
659 values<0, false, false, true>(
660 get_global_shape_values<Number>(mf_object_id), u, u);
661 __syncthreads();
662
663 break;
664 }
665 case 3:
666 {
667 gradients<2, false, true, false>(
668 get_global_co_shape_gradients<Number>(mf_object_id),
669 grad_u[2],
670 u);
671 __syncthreads();
672 gradients<1, false, true, false>(
673 get_global_co_shape_gradients<Number>(mf_object_id),
674 grad_u[1],
675 u);
676 __syncthreads();
677 gradients<0, false, true, false>(
678 get_global_co_shape_gradients<Number>(mf_object_id),
679 grad_u[0],
680 u);
681 __syncthreads();
682
683 values<2, false, false, true>(
684 get_global_shape_values<Number>(mf_object_id), u, u);
685 __syncthreads();
686 values<1, false, false, true>(
687 get_global_shape_values<Number>(mf_object_id), u, u);
688 __syncthreads();
689 values<0, false, false, true>(
690 get_global_shape_values<Number>(mf_object_id), u, u);
691 __syncthreads();
692
693 break;
694 }
695 default:
696 {
697 // Do nothing. We should throw but we can't from a __device__
698 // function.
699 }
700 }
701 }
702 } // namespace internal
703} // namespace CUDAWrappers
704
705DEAL_II_NAMESPACE_CLOSE
706
707#endif
708
709#endif
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
CUDAWrappers::internal::EvaluatorVariant
EvaluatorVariant
Definition: cuda_tensor_product_kernels.h:43
CUDAWrappers::internal::evaluate_general
@ evaluate_general
Definition: cuda_tensor_product_kernels.h:44
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
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
CUDAWrappers::internal::EvaluatorTensorProduct< evaluate_general, dim, fe_degree, n_q_points_1d, Number >::mf_object_id
const int mf_object_id
Definition: cuda_tensor_product_kernels.h:153
Generated by doxygen 1.9.4