Reference documentation for deal.II version Git 0994510114 2020-11-30 18:32:56 -0700
\(\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 
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) ?
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
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
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
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
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
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
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 
706 
707 #endif
708 
709 #endif
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:461
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371