Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
cuda_tensor_product_kernels.h
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
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>
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__
86 
91  template <int direction, bool dof_to_quad, bool add, bool in_place>
92  __device__ void
93  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(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  template <bool add>
134  __device__ void
135  integrate_gradient(Number *u, Number *grad_u[dim]);
136  };
137 
138 
139 
140  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
141  __device__
142  EvaluatorTensorProduct<evaluate_general,
143  dim,
144  fe_degree,
145  n_q_points_1d,
146  Number>::EvaluatorTensorProduct()
147  {}
148 
149 
150 
151  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
152  template <int direction, bool dof_to_quad, bool add, bool in_place>
153  __device__ void
154  EvaluatorTensorProduct<evaluate_general,
155  dim,
156  fe_degree,
157  n_q_points_1d,
158  Number>::values(const Number *in, Number *out) const
159  {
160  apply<direction, dof_to_quad, add, in_place>(global_shape_values,
161  in,
162  out);
163  }
164 
165 
166 
167  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
168  template <int direction, bool dof_to_quad, bool add, bool in_place>
169  __device__ void
170  EvaluatorTensorProduct<evaluate_general,
171  dim,
172  fe_degree,
173  n_q_points_1d,
174  Number>::gradients(const Number *in,
175  Number * out) const
176  {
177  apply<direction, dof_to_quad, add, in_place>(global_shape_gradients,
178  in,
179  out);
180  }
181 
182 
183 
184  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
185  template <int direction, bool dof_to_quad, bool add, bool in_place>
186  __device__ void
187  EvaluatorTensorProduct<evaluate_general,
188  dim,
189  fe_degree,
190  n_q_points_1d,
191  Number>::apply(Number shape_data[],
192  const Number *in,
193  Number * out) const
194  {
195  const unsigned int i = (dim == 1) ? 0 : threadIdx.x % n_q_points_1d;
196  const unsigned int j = (dim == 3) ? threadIdx.y : 0;
197  const unsigned int q = (dim == 1) ?
198  (threadIdx.x % n_q_points_1d) :
199  (dim == 2) ? threadIdx.y : threadIdx.z;
200 
201  // This loop simply multiply the shape function at the quadrature point by
202  // the value finite element coefficient.
203  Number t = 0;
204  for (int k = 0; k < n_q_points_1d; ++k)
205  {
206  const unsigned int shape_idx =
207  dof_to_quad ? (q + k * n_q_points_1d) : (k + q * n_q_points_1d);
208  const unsigned int source_idx =
209  (direction == 0) ?
210  (k + n_q_points_1d * (i + n_q_points_1d * j)) :
211  (direction == 1) ? (i + n_q_points_1d * (k + n_q_points_1d * j)) :
212  (i + n_q_points_1d * (j + n_q_points_1d * k));
213  t += shape_data[shape_idx] *
214  (in_place ? out[source_idx] : in[source_idx]);
215  }
216 
217  if (in_place)
218  __syncthreads();
219 
220  const unsigned int destination_idx =
221  (direction == 0) ?
222  (q + n_q_points_1d * (i + n_q_points_1d * j)) :
223  (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) :
224  (i + n_q_points_1d * (j + n_q_points_1d * q));
225 
226  if (add)
227  out[destination_idx] += t;
228  else
229  out[destination_idx] = t;
230  }
231 
232 
233 
234  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
235  inline __device__ void
236  EvaluatorTensorProduct<evaluate_general,
237  dim,
238  fe_degree,
239  n_q_points_1d,
240  Number>::value_at_quad_pts(Number *u)
241  {
242  switch (dim)
243  {
244  case 1:
245  {
246  values<0, true, false, true>(u, u);
247 
248  break;
249  }
250  case 2:
251  {
252  values<0, true, false, true>(u, u);
253  __syncthreads();
254  values<1, true, false, true>(u, u);
255 
256  break;
257  }
258  case 3:
259  {
260  values<0, true, false, true>(u, u);
261  __syncthreads();
262  values<1, true, false, true>(u, u);
263  __syncthreads();
264  values<2, true, false, true>(u, u);
265 
266  break;
267  }
268  default:
269  {
270  // Do nothing. We should throw but we can't from a __device__
271  // function.
272  }
273  }
274  }
275 
276 
277 
278  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
279  inline __device__ void
280  EvaluatorTensorProduct<evaluate_general,
281  dim,
282  fe_degree,
283  n_q_points_1d,
284  Number>::integrate_value(Number *u)
285  {
286  switch (dim)
287  {
288  case 1:
289  {
290  values<0, false, false, true>(u, u);
291 
292  break;
293  }
294  case 2:
295  {
296  values<0, false, false, true>(u, u);
297  __syncthreads();
298  values<1, false, false, true>(u, u);
299 
300  break;
301  }
302  case 3:
303  {
304  values<0, false, false, true>(u, u);
305  __syncthreads();
306  values<1, false, false, true>(u, u);
307  __syncthreads();
308  values<2, false, false, true>(u, u);
309 
310  break;
311  }
312  default:
313  {
314  // Do nothing. We should throw but we can't from a __device__
315  // function.
316  }
317  }
318  }
319 
320 
321 
322  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
323  inline __device__ void
324  EvaluatorTensorProduct<evaluate_general,
325  dim,
326  fe_degree,
327  n_q_points_1d,
328  Number>::gradient_at_quad_pts(const Number *const u,
329  Number *grad_u[dim])
330  {
331  switch (dim)
332  {
333  case 1:
334  {
335  gradients<0, true, false, false>(u, grad_u[0]);
336 
337  break;
338  }
339  case 2:
340  {
341  gradients<0, true, false, false>(u, grad_u[0]);
342  values<0, true, false, false>(u, grad_u[1]);
343 
344  __syncthreads();
345 
346  values<1, true, false, true>(grad_u[0], grad_u[0]);
347  gradients<1, true, false, true>(grad_u[1], grad_u[1]);
348 
349  break;
350  }
351  case 3:
352  {
353  gradients<0, true, false, false>(u, grad_u[0]);
354  values<0, true, false, false>(u, grad_u[1]);
355  values<0, true, false, false>(u, grad_u[2]);
356 
357  __syncthreads();
358 
359  values<1, true, false, true>(grad_u[0], grad_u[0]);
360  gradients<1, true, false, true>(grad_u[1], grad_u[1]);
361  values<1, true, false, true>(grad_u[2], grad_u[2]);
362 
363  __syncthreads();
364 
365  values<2, true, false, true>(grad_u[0], grad_u[0]);
366  values<2, true, false, true>(grad_u[1], grad_u[1]);
367  gradients<2, true, false, true>(grad_u[2], grad_u[2]);
368 
369  break;
370  }
371  default:
372  {
373  // Do nothing. We should throw but we can't from a __device__
374  // function.
375  }
376  }
377  }
378 
379 
380 
381  template <int dim, int fe_degree, int n_q_points_1d, typename Number>
382  template <bool add>
383  inline __device__ void
384  EvaluatorTensorProduct<evaluate_general,
385  dim,
386  fe_degree,
387  n_q_points_1d,
388  Number>::integrate_gradient(Number *u,
389  Number *grad_u[dim])
390  {
391  switch (dim)
392  {
393  case 1:
394  {
395  gradients<0, false, add, false>(grad_u[dim], u);
396 
397  break;
398  }
399  case 2:
400  {
401  gradients<0, false, false, true>(grad_u[0], grad_u[0]);
402  values<0, false, false, true>(grad_u[1], grad_u[1]);
403 
404  __syncthreads();
405 
406  values<1, false, add, false>(grad_u[0], u);
407  __syncthreads();
408  gradients<1, false, true, false>(grad_u[1], u);
409 
410  break;
411  }
412  case 3:
413  {
414  gradients<0, false, false, true>(grad_u[0], grad_u[0]);
415  values<0, false, false, true>(grad_u[1], grad_u[1]);
416  values<0, false, false, true>(grad_u[2], grad_u[2]);
417 
418  __syncthreads();
419 
420  values<1, false, false, true>(grad_u[0], grad_u[0]);
421  gradients<1, false, false, true>(grad_u[1], grad_u[1]);
422  values<1, false, false, true>(grad_u[2], grad_u[2]);
423 
424  __syncthreads();
425 
426  values<2, false, add, false>(grad_u[0], u);
427  __syncthreads();
428  values<2, false, true, false>(grad_u[1], u);
429  __syncthreads();
430  gradients<2, false, true, false>(grad_u[2], u);
431 
432  break;
433  }
434  default:
435  {
436  // Do nothing. We should throw but we can't from a __device__
437  // function.
438  }
439  }
440  }
441  } // namespace internal
442 } // namespace CUDAWrappers
443 
444 DEAL_II_NAMESPACE_CLOSE
445 
446 #endif
447 
448 #endif
constexpr unsigned int pow(const unsigned int base, const int iexp)
Definition: utilities.h:428