Reference documentation for deal.II version GIT c14369f203 2023-10-01 07:40:02+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_fe_evaluation.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2016 - 2023 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 #ifndef dealii_cuda_fe_evaluation_h
17 #define dealii_cuda_fe_evaluation_h
18 
19 #include <deal.II/base/config.h>
20 
22 #include <deal.II/base/tensor.h>
23 #include <deal.II/base/utilities.h>
24 
27 
30 #include <deal.II/matrix_free/cuda_matrix_free.templates.h>
32 
33 #include <Kokkos_Core.hpp>
34 
36 
40 namespace CUDAWrappers
41 {
67  template <int dim,
68  int fe_degree,
69  int n_q_points_1d = fe_degree + 1,
70  int n_components_ = 1,
71  typename Number = double>
73  {
74  public:
78  using value_type = Number;
79 
84 
89 
93  static constexpr unsigned int dimension = dim;
94 
98  static constexpr unsigned int n_components = n_components_;
99 
103  static constexpr unsigned int n_q_points =
104  Utilities::pow(n_q_points_1d, dim);
105 
109  static constexpr unsigned int tensor_dofs_per_cell =
110  Utilities::pow(fe_degree + 1, dim);
111 
117 
127  read_dof_values(const Number *src);
128 
136  distribute_local_to_global(Number *dst) const;
137 
146  evaluate(const bool evaluate_val, const bool evaluate_grad);
147 
156  integrate(const bool integrate_val, const bool integrate_grad);
157 
163  get_value(int q_point) const;
164 
170  get_dof_value(int q_point) const;
171 
177  submit_value(const value_type &val_in, int q_point);
178 
184  submit_dof_value(const value_type &val_in, int q_point);
185 
191  get_gradient(int q_point) const;
192 
198  submit_gradient(const gradient_type &grad_in, int q_point);
199 
200  // clang-format off
211  // clang-format on
212  template <typename Functor>
214  apply_for_each_quad_point(const Functor &func);
215 
216  private:
217  const data_type *data;
219  int cell_id;
220  };
221 
222 
223 
224  template <int dim,
225  int fe_degree,
226  int n_q_points_1d,
227  int n_components_,
228  typename Number>
231  FEEvaluation(const data_type *data, SharedData<dim, Number> *shdata)
232  : data(data)
233  , shared_data(shdata)
234  , cell_id(shared_data->team_member.league_rank())
235  {}
236 
237 
238 
239  template <int dim,
240  int fe_degree,
241  int n_q_points_1d,
242  int n_components_,
243  typename Number>
246  read_dof_values(const Number *src)
247  {
248  static_assert(n_components_ == 1, "This function only supports FE with one \
249  components");
250  // Populate the scratch memory
251  Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
252  n_q_points),
253  [&](const int &i) {
254  shared_data->values(i) =
255  src[data->local_to_global(cell_id, i)];
256  });
257  shared_data->team_member.team_barrier();
258 
259  internal::resolve_hanging_nodes<dim, fe_degree, false>(
260  shared_data->team_member,
261  data->constraint_weights,
262  data->constraint_mask(cell_id),
263  shared_data->values);
264  }
265 
266 
267 
268  template <int dim,
269  int fe_degree,
270  int n_q_points_1d,
271  int n_components_,
272  typename Number>
275  distribute_local_to_global(Number *dst) const
276  {
277  static_assert(n_components_ == 1, "This function only supports FE with one \
278  components");
279 
280  internal::resolve_hanging_nodes<dim, fe_degree, true>(
281  shared_data->team_member,
282  data->constraint_weights,
283  data->constraint_mask(cell_id),
284  shared_data->values);
285 
286  if (data->use_coloring)
287  {
288  Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
289  n_q_points),
290  [&](const int &i) {
291  dst[data->local_to_global(cell_id, i)] +=
292  shared_data->values(i);
293  });
294  }
295  else
296  {
298  Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
299  [&](const int &i) {
300  Kokkos::atomic_add(&dst[data->local_to_global(cell_id, i)],
301  shared_data->values(i));
302  });
303  }
304  }
305 
306 
307 
308  template <int dim,
309  int fe_degree,
310  int n_q_points_1d,
311  int n_components_,
312  typename Number>
315  const bool evaluate_val,
316  const bool evaluate_grad)
317  {
318  // First evaluate the gradients because it requires values that will be
319  // changed if evaluate_val is true
322  dim,
323  fe_degree,
324  n_q_points_1d,
325  Number>
326  evaluator_tensor_product(shared_data->team_member,
327  data->shape_values,
328  data->shape_gradients,
329  data->co_shape_gradients);
330  if (evaluate_val == true && evaluate_grad == true)
331  {
332  evaluator_tensor_product.value_and_gradient_at_quad_pts(
333  shared_data->values, shared_data->gradients);
334  shared_data->team_member.team_barrier();
335  }
336  else if (evaluate_grad == true)
337  {
338  evaluator_tensor_product.gradient_at_quad_pts(shared_data->values,
339  shared_data->gradients);
340  shared_data->team_member.team_barrier();
341  }
342  else if (evaluate_val == true)
343  {
344  evaluator_tensor_product.value_at_quad_pts(shared_data->values);
345  shared_data->team_member.team_barrier();
346  }
347  }
348 
349 
350 
351  template <int dim,
352  int fe_degree,
353  int n_q_points_1d,
354  int n_components_,
355  typename Number>
358  const bool integrate_val,
359  const bool integrate_grad)
360  {
363  dim,
364  fe_degree,
365  n_q_points_1d,
366  Number>
367  evaluator_tensor_product(shared_data->team_member,
368  data->shape_values,
369  data->shape_gradients,
370  data->co_shape_gradients);
371  if (integrate_val == true && integrate_grad == true)
372  {
373  evaluator_tensor_product.integrate_value_and_gradient(
374  shared_data->values, shared_data->gradients);
375  }
376  else if (integrate_val == true)
377  {
378  evaluator_tensor_product.integrate_value(shared_data->values);
379  shared_data->team_member.team_barrier();
380  }
381  else if (integrate_grad == true)
382  {
383  evaluator_tensor_product.template integrate_gradient<false>(
384  shared_data->values, shared_data->gradients);
385  shared_data->team_member.team_barrier();
386  }
387  }
388 
389 
390 
391  template <int dim,
392  int fe_degree,
393  int n_q_points_1d,
394  int n_components_,
395  typename Number>
396  DEAL_II_HOST_DEVICE typename FEEvaluation<dim,
397  fe_degree,
398  n_q_points_1d,
399  n_components_,
400  Number>::value_type
402  int q_point) const
403  {
404  return shared_data->values(q_point);
405  }
406 
407 
408 
409  template <int dim,
410  int fe_degree,
411  int n_q_points_1d,
412  int n_components_,
413  typename Number>
414  DEAL_II_HOST_DEVICE typename FEEvaluation<dim,
415  fe_degree,
416  n_q_points_1d,
417  n_components_,
418  Number>::value_type
420  get_dof_value(int q_point) const
421  {
422  return shared_data->values(q_point);
423  }
424 
425 
426 
427  template <int dim,
428  int fe_degree,
429  int n_q_points_1d,
430  int n_components_,
431  typename Number>
434  submit_value(const value_type &val_in, int q_point)
435  {
436  shared_data->values(q_point) = val_in * data->JxW(cell_id, q_point);
437  }
438 
439 
440 
441  template <int dim,
442  int fe_degree,
443  int n_q_points_1d,
444  int n_components_,
445  typename Number>
448  submit_dof_value(const value_type &val_in, int q_point)
449  {
450  shared_data->values(q_point) = val_in;
451  }
452 
453 
454 
455  template <int dim,
456  int fe_degree,
457  int n_q_points_1d,
458  int n_components_,
459  typename Number>
460  DEAL_II_HOST_DEVICE typename FEEvaluation<dim,
461  fe_degree,
462  n_q_points_1d,
463  n_components_,
464  Number>::gradient_type
466  get_gradient(int q_point) const
467  {
468  static_assert(n_components_ == 1, "This function only supports FE with one \
469  components");
470 
471  gradient_type grad;
472  for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
473  {
474  Number tmp = 0.;
475  for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
476  tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) *
477  shared_data->gradients(q_point, d_2);
478  grad[d_1] = tmp;
479  }
480 
481  return grad;
482  }
483 
484 
485 
486  template <int dim,
487  int fe_degree,
488  int n_q_points_1d,
489  int n_components_,
490  typename Number>
493  submit_gradient(const gradient_type &grad_in, int q_point)
494  {
495  for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
496  {
497  Number tmp = 0.;
498  for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
499  tmp += data->inv_jacobian(cell_id, q_point, d_1, d_2) * grad_in[d_2];
500  shared_data->gradients(q_point, d_1) =
501  tmp * data->JxW(cell_id, q_point);
502  }
503  }
504 
505 
506 
507  template <int dim,
508  int fe_degree,
509  int n_q_points_1d,
510  int n_components_,
511  typename Number>
512  template <typename Functor>
515  apply_for_each_quad_point(const Functor &func)
516  {
517  Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
518  n_q_points),
519  [&](const int &i) { func(this, i); });
520  shared_data->team_member.team_barrier();
521  }
522 
523 
524 
525 #ifndef DOXYGEN
526  template <int dim,
527  int fe_degree,
528  int n_q_points_1d,
529  int n_components_,
530  typename Number>
531  constexpr unsigned int
533  n_q_points;
534 #endif
535 } // namespace CUDAWrappers
536 
538 
539 #endif
static constexpr unsigned int tensor_dofs_per_cell
void integrate(const bool integrate_val, const bool integrate_grad)
FEEvaluation(const data_type *data, SharedData< dim, Number > *shdata)
value_type get_value(int q_point) const
value_type get_dof_value(int q_point) const
SharedData< dim, Number > * shared_data
void distribute_local_to_global(Number *dst) const
static constexpr unsigned int dimension
void evaluate(const bool evaluate_val, const bool evaluate_grad)
static constexpr unsigned int n_q_points
gradient_type get_gradient(int q_point) const
static constexpr unsigned int n_components
typename MatrixFree< dim, Number >::Data data_type
void read_dof_values(const Number *src)
void submit_dof_value(const value_type &val_in, int q_point)
void submit_value(const value_type &val_in, int q_point)
void submit_gradient(const gradient_type &grad_in, int q_point)
void apply_for_each_quad_point(const Functor &func)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:447
void parallel_for(Iterator x_begin, Iterator x_end, const Functor &functor, const unsigned int grainsize)
Definition: parallel.h:84
#define DEAL_II_HOST_DEVICE
Definition: numbers.h:35