Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
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
23
24#include <deal.II/matrix_free/cuda_matrix_free.templates.h>
25
26#ifdef DEAL_II_COMPILER_CUDA_AWARE
27
29
30
31namespace 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:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:462