Reference documentation for deal.II version 9.2.0
\(\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 - 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 #ifndef dealii_cuda_fe_evaluation_h
17 #define dealii_cuda_fe_evaluation_h
18 
19 #include <deal.II/base/config.h>
20 
21 #ifdef DEAL_II_COMPILER_CUDA_AWARE
22 
23 # include <deal.II/base/tensor.h>
24 # include <deal.II/base/utilities.h>
25 
26 # include <deal.II/lac/cuda_atomic.h>
27 # include <deal.II/lac/cuda_vector.h>
28 
29 # include <deal.II/matrix_free/cuda_hanging_nodes_internal.h>
30 # include <deal.II/matrix_free/cuda_matrix_free.h>
31 # include <deal.II/matrix_free/cuda_matrix_free.templates.h>
32 # include <deal.II/matrix_free/cuda_tensor_product_kernels.h>
33 
34 DEAL_II_NAMESPACE_OPEN
35 
39 namespace CUDAWrappers
40 {
41  namespace internal
42  {
47  template <int dim, int n_points_1d>
48  __device__ inline unsigned int
49  compute_index()
50  {
51  return (dim == 1 ?
52  threadIdx.x % n_points_1d :
53  dim == 2 ?
54  threadIdx.x % n_points_1d + n_points_1d * threadIdx.y :
55  threadIdx.x % n_points_1d +
56  n_points_1d * (threadIdx.y + n_points_1d * threadIdx.z));
57  }
58  } // namespace internal
59 
87  template <int dim,
88  int fe_degree,
89  int n_q_points_1d = fe_degree + 1,
90  int n_components_ = 1,
91  typename Number = double>
92  class FEEvaluation
93  {
94  public:
98  using value_type = Number;
99 
103  using gradient_type = Tensor<1, dim, Number>;
104 
108  using data_type = typename MatrixFree<dim, Number>::Data;
109 
113  static constexpr unsigned int dimension = dim;
114 
118  static constexpr unsigned int n_components = n_components_;
119 
123  static constexpr unsigned int n_q_points =
124  Utilities::pow(n_q_points_1d, dim);
125 
129  static constexpr unsigned int tensor_dofs_per_cell =
130  Utilities::pow(fe_degree + 1, dim);
131 
135  __device__
136  FEEvaluation(const unsigned int cell_id,
137  const data_type * data,
138  SharedData<dim, Number> *shdata);
139 
148  __device__ void
149  read_dof_values(const Number *src);
150 
157  __device__ void
158  distribute_local_to_global(Number *dst) const;
159 
167  __device__ void
168  evaluate(const bool evaluate_val, const bool evaluate_grad);
169 
177  __device__ void
178  integrate(const bool integrate_val, const bool integrate_grad);
179 
186  DEAL_II_DEPRECATED __device__ value_type
187  get_value(const unsigned int q_point) const;
188 
193  __device__ value_type
194  get_value() const;
195 
202  DEAL_II_DEPRECATED __device__ value_type
203  get_dof_value(const unsigned int dof) const;
204 
209  __device__ value_type
210  get_dof_value() const;
211 
220  DEAL_II_DEPRECATED __device__ void
221  submit_value(const value_type &val_in, const unsigned int q_point);
222 
227  __device__ void
228  submit_value(const value_type &val_in);
229 
238  DEAL_II_DEPRECATED __device__ void
239  submit_dof_value(const value_type &val_in, const unsigned int dof);
240 
245  __device__ void
246  submit_dof_value(const value_type &val_in);
247 
254  DEAL_II_DEPRECATED __device__ gradient_type
255  get_gradient(const unsigned int q_point) const;
256 
261  __device__ gradient_type
262  get_gradient() const;
263 
270  DEAL_II_DEPRECATED __device__ void
271  submit_gradient(const gradient_type &grad_in, const unsigned int q_point);
272 
273 
278  __device__ void
279  submit_gradient(const gradient_type &grad_in);
280 
281  // clang-format off
294  // clang-format on
295  template <typename Functor>
296  DEAL_II_DEPRECATED __device__ void
297  apply_quad_point_operations(const Functor &func);
298 
299  // clang-format off
310  // clang-format on
311  template <typename Functor>
312  __device__ void
313  apply_for_each_quad_point(const Functor &func);
314 
315  private:
316  types::global_dof_index *local_to_global;
317  unsigned int n_cells;
318  unsigned int padding_length;
319 
320  const unsigned int constraint_mask;
321 
322  const bool use_coloring;
323 
324  Number *inv_jac;
325  Number *JxW;
326 
327  // Internal buffer
328  Number *values;
329  Number *gradients[dim];
330  };
331 
332 
333 
334  template <int dim,
335  int fe_degree,
336  int n_q_points_1d,
337  int n_components_,
338  typename Number>
339  __device__
340  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
341  FEEvaluation(const unsigned int cell_id,
342  const data_type * data,
343  SharedData<dim, Number> *shdata)
344  : n_cells(data->n_cells)
345  , padding_length(data->padding_length)
346  , constraint_mask(data->constraint_mask[cell_id])
347  , use_coloring(data->use_coloring)
348  , values(shdata->values)
349  {
350  local_to_global = data->local_to_global + padding_length * cell_id;
351  inv_jac = data->inv_jacobian + padding_length * cell_id;
352  JxW = data->JxW + padding_length * cell_id;
353 
354  for (unsigned int i = 0; i < dim; ++i)
355  gradients[i] = shdata->gradients[i];
356  }
357 
358 
359 
360  template <int dim,
361  int fe_degree,
362  int n_q_points_1d,
363  int n_components_,
364  typename Number>
365  __device__ void
366  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
367  read_dof_values(const Number *src)
368  {
369  static_assert(n_components_ == 1, "This function only supports FE with one \
370  components");
371  const unsigned int idx = internal::compute_index<dim, n_q_points_1d>();
372 
373  const types::global_dof_index src_idx = local_to_global[idx];
374  // Use the read-only data cache.
375  values[idx] = __ldg(&src[src_idx]);
376 
377  __syncthreads();
378 
379  internal::resolve_hanging_nodes<dim, fe_degree, false>(constraint_mask,
380  values);
381  }
382 
383 
384 
385  template <int dim,
386  int fe_degree,
387  int n_q_points_1d,
388  int n_components_,
389  typename Number>
390  __device__ void
391  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
392  distribute_local_to_global(Number *dst) const
393  {
394  static_assert(n_components_ == 1, "This function only supports FE with one \
395  components");
396  internal::resolve_hanging_nodes<dim, fe_degree, true>(constraint_mask,
397  values);
398 
399  const unsigned int idx = internal::compute_index<dim, n_q_points_1d>();
400 
401  const types::global_dof_index destination_idx = local_to_global[idx];
402 
403  if (use_coloring)
404  dst[destination_idx] += values[idx];
405  else
406  LinearAlgebra::CUDAWrappers::atomicAdd_wrapper(&dst[destination_idx],
407  values[idx]);
408  }
409 
410 
411 
412  template <int dim,
413  int fe_degree,
414  int n_q_points_1d,
415  int n_components_,
416  typename Number>
417  __device__ void
418  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::evaluate(
419  const bool evaluate_val,
420  const bool evaluate_grad)
421  {
422  // First evaluate the gradients because it requires values that will be
423  // changed if evaluate_val is true
424  internal::EvaluatorTensorProduct<
425  internal::EvaluatorVariant::evaluate_general,
426  dim,
427  fe_degree,
428  n_q_points_1d,
429  Number>
430  evaluator_tensor_product;
431  if (evaluate_val == true && evaluate_grad == true)
432  {
433  evaluator_tensor_product.value_and_gradient_at_quad_pts(values,
434  gradients);
435  __syncthreads();
436  }
437  else if (evaluate_grad == true)
438  {
439  evaluator_tensor_product.gradient_at_quad_pts(values, gradients);
440  __syncthreads();
441  }
442  else if (evaluate_val == true)
443  {
444  evaluator_tensor_product.value_at_quad_pts(values);
445  __syncthreads();
446  }
447  }
448 
449 
450 
451  template <int dim,
452  int fe_degree,
453  int n_q_points_1d,
454  int n_components_,
455  typename Number>
456  __device__ void
457  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::integrate(
458  const bool integrate_val,
459  const bool integrate_grad)
460  {
461  internal::EvaluatorTensorProduct<
462  internal::EvaluatorVariant::evaluate_general,
463  dim,
464  fe_degree,
465  n_q_points_1d,
466  Number>
467  evaluator_tensor_product;
468  if (integrate_val == true && integrate_grad == true)
469  {
470  evaluator_tensor_product.integrate_value_and_gradient(values,
471  gradients);
472  }
473  else if (integrate_val == true)
474  {
475  evaluator_tensor_product.integrate_value(values);
476  __syncthreads();
477  }
478  else if (integrate_grad == true)
479  {
480  evaluator_tensor_product.integrate_gradient<false>(values, gradients);
481  __syncthreads();
482  }
483  }
484 
485 
486 
487  template <int dim,
488  int fe_degree,
489  int n_q_points_1d,
490  int n_components_,
491  typename Number>
492  __device__ typename FEEvaluation<dim,
493  fe_degree,
494  n_q_points_1d,
495  n_components_,
496  Number>::value_type
497  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::get_value(
498  const unsigned int q_point) const
499  {
500  return values[q_point];
501  }
502 
503 
504 
505  template <int dim,
506  int fe_degree,
507  int n_q_points_1d,
508  int n_components_,
509  typename Number>
510  __device__ typename FEEvaluation<dim,
511  fe_degree,
512  n_q_points_1d,
513  n_components_,
514  Number>::value_type
515  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
516  get_value() const
517  {
518  const unsigned int q_point = internal::compute_index<dim, n_q_points_1d>();
519  return values[q_point];
520  }
521 
522 
523 
524  template <int dim,
525  int fe_degree,
526  int n_q_points_1d,
527  int n_components_,
528  typename Number>
529  __device__ typename FEEvaluation<dim,
530  fe_degree,
531  n_q_points_1d,
532  n_components_,
533  Number>::value_type
534  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
535  get_dof_value(const unsigned int dof) const
536  {
537  return values[dof];
538  }
539 
540 
541 
542  template <int dim,
543  int fe_degree,
544  int n_q_points_1d,
545  int n_components_,
546  typename Number>
547  __device__ typename FEEvaluation<dim,
548  fe_degree,
549  n_q_points_1d,
550  n_components_,
551  Number>::value_type
552  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
553  get_dof_value() const
554  {
555  const unsigned int dof = internal::compute_index<dim, fe_degree + 1>();
556  return values[dof];
557  }
558 
559 
560 
561  template <int dim,
562  int fe_degree,
563  int n_q_points_1d,
564  int n_components_,
565  typename Number>
566  __device__ void
567  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
568  submit_value(const value_type &val_in, const unsigned int q_point)
569  {
570  values[q_point] = val_in * JxW[q_point];
571  }
572 
573 
574 
575  template <int dim,
576  int fe_degree,
577  int n_q_points_1d,
578  int n_components_,
579  typename Number>
580  __device__ void
581  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
582  submit_value(const value_type &val_in)
583  {
584  const unsigned int q_point = internal::compute_index<dim, n_q_points_1d>();
585  values[q_point] = val_in * JxW[q_point];
586  }
587 
588 
589 
590  template <int dim,
591  int fe_degree,
592  int n_q_points_1d,
593  int n_components_,
594  typename Number>
595  __device__ void
596  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
597  submit_dof_value(const value_type &val_in, const unsigned int dof)
598  {
599  values[dof] = val_in;
600  }
601 
602 
603 
604  template <int dim,
605  int fe_degree,
606  int n_q_points_1d,
607  int n_components_,
608  typename Number>
609  __device__ void
610  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
611  submit_dof_value(const value_type &val_in)
612  {
613  const unsigned int dof = internal::compute_index<dim, fe_degree + 1>();
614  values[dof] = val_in;
615  }
616 
617 
618 
619  template <int dim,
620  int fe_degree,
621  int n_q_points_1d,
622  int n_components_,
623  typename Number>
624  __device__ typename FEEvaluation<dim,
625  fe_degree,
626  n_q_points_1d,
627  n_components_,
628  Number>::gradient_type
629  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
630  get_gradient(const unsigned int q_point) const
631  {
632  static_assert(n_components_ == 1, "This function only supports FE with one \
633  components");
634  // TODO optimize if the mesh is uniform
635  const Number *inv_jacobian = &inv_jac[q_point];
636  gradient_type grad;
637  for (int d_1 = 0; d_1 < dim; ++d_1)
638  {
639  Number tmp = 0.;
640  for (int d_2 = 0; d_2 < dim; ++d_2)
641  tmp += inv_jacobian[padding_length * n_cells * (dim * d_2 + d_1)] *
642  gradients[d_2][q_point];
643  grad[d_1] = tmp;
644  }
645 
646  return grad;
647  }
648 
649 
650 
651  template <int dim,
652  int fe_degree,
653  int n_q_points_1d,
654  int n_components_,
655  typename Number>
656  __device__ typename FEEvaluation<dim,
657  fe_degree,
658  n_q_points_1d,
659  n_components_,
660  Number>::gradient_type
661  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
662  get_gradient() const
663  {
664  static_assert(n_components_ == 1, "This function only supports FE with one \
665  components");
666 
667  // TODO optimize if the mesh is uniform
668  const unsigned int q_point = internal::compute_index<dim, n_q_points_1d>();
669  const Number * inv_jacobian = &inv_jac[q_point];
670  gradient_type grad;
671  for (int d_1 = 0; d_1 < dim; ++d_1)
672  {
673  Number tmp = 0.;
674  for (int d_2 = 0; d_2 < dim; ++d_2)
675  tmp += inv_jacobian[padding_length * n_cells * (dim * d_2 + d_1)] *
676  gradients[d_2][q_point];
677  grad[d_1] = tmp;
678  }
679 
680  return grad;
681  }
682 
683 
684 
685  template <int dim,
686  int fe_degree,
687  int n_q_points_1d,
688  int n_components_,
689  typename Number>
690  __device__ void
691  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
692  submit_gradient(const gradient_type &grad_in, const unsigned int q_point)
693  {
694  // TODO optimize if the mesh is uniform
695  const Number *inv_jacobian = &inv_jac[q_point];
696  for (int d_1 = 0; d_1 < dim; ++d_1)
697  {
698  Number tmp = 0.;
699  for (int d_2 = 0; d_2 < dim; ++d_2)
700  tmp += inv_jacobian[n_cells * padding_length * (dim * d_1 + d_2)] *
701  grad_in[d_2];
702  gradients[d_1][q_point] = tmp * JxW[q_point];
703  }
704  }
705 
706 
707 
708  template <int dim,
709  int fe_degree,
710  int n_q_points_1d,
711  int n_components_,
712  typename Number>
713  __device__ void
714  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
715  submit_gradient(const gradient_type &grad_in)
716  {
717  // TODO optimize if the mesh is uniform
718  const unsigned int q_point = internal::compute_index<dim, n_q_points_1d>();
719  const Number * inv_jacobian = &inv_jac[q_point];
720  for (int d_1 = 0; d_1 < dim; ++d_1)
721  {
722  Number tmp = 0.;
723  for (int d_2 = 0; d_2 < dim; ++d_2)
724  tmp += inv_jacobian[n_cells * padding_length * (dim * d_1 + d_2)] *
725  grad_in[d_2];
726  gradients[d_1][q_point] = tmp * JxW[q_point];
727  }
728  }
729 
730 
731 
732  template <int dim,
733  int fe_degree,
734  int n_q_points_1d,
735  int n_components_,
736  typename Number>
737  template <typename Functor>
738  __device__ void
739  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
740  apply_quad_point_operations(const Functor &func)
741  {
742  func(this, internal::compute_index<dim, n_q_points_1d>());
743 
744  __syncthreads();
745  }
746 
747 
748 
749  template <int dim,
750  int fe_degree,
751  int n_q_points_1d,
752  int n_components_,
753  typename Number>
754  template <typename Functor>
755  __device__ void
756  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
757  apply_for_each_quad_point(const Functor &func)
758  {
759  func(this);
760 
761  __syncthreads();
762  }
763 } // namespace CUDAWrappers
764 
765 DEAL_II_NAMESPACE_CLOSE
766 
767 #endif
768 
769 #endif
CUDAWrappers::FEEvaluation::dimension
static constexpr unsigned int dimension
Definition: cuda_fe_evaluation.h:113
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
CUDAWrappers::FEEvaluation::local_to_global
types::global_dof_index * local_to_global
Definition: cuda_fe_evaluation.h:316
CUDAWrappers::internal::evaluate_general
@ evaluate_general
Definition: cuda_tensor_product_kernels.h:44
cuda_matrix_free.h
CUDAWrappers::FEEvaluation::data_type
typename MatrixFree< dim, Number >::Data data_type
Definition: cuda_fe_evaluation.h:108
CUDAWrappers::FEEvaluation::submit_dof_value
void submit_dof_value(const value_type &val_in, const unsigned int dof)
Definition: cuda_fe_evaluation.h:597
cuda_vector.h
CUDAWrappers::FEEvaluation::FEEvaluation
FEEvaluation(const unsigned int cell_id, const data_type *data, SharedData< dim, Number > *shdata)
Definition: cuda_fe_evaluation.h:341
CUDAWrappers::FEEvaluation::n_components
static constexpr unsigned int n_components
Definition: cuda_fe_evaluation.h:118
CUDAWrappers::SharedData::gradients
Number * gradients[dim]
Definition: cuda_matrix_free.h:614
utilities.h
CUDAWrappers::FEEvaluation::get_dof_value
value_type get_dof_value() const
Definition: cuda_fe_evaluation.h:553
CUDAWrappers::FEEvaluation::padding_length
unsigned int padding_length
Definition: cuda_fe_evaluation.h:318
LinearAlgebra::CUDAWrappers::atomicAdd_wrapper
float atomicAdd_wrapper(float *address, float val)
Definition: cuda_atomic.h:35
CUDAWrappers::SharedData
Definition: cuda_matrix_free.h:591
CUDAWrappers::FEEvaluation::read_dof_values
void read_dof_values(const Number *src)
Definition: cuda_fe_evaluation.h:367
CUDAWrappers::FEEvaluation::apply_quad_point_operations
void apply_quad_point_operations(const Functor &func)
Definition: cuda_fe_evaluation.h:740
CUDAWrappers::FEEvaluation::value_type
Number value_type
Definition: cuda_fe_evaluation.h:98
CUDAWrappers::internal::EvaluatorTensorProduct
Definition: cuda_tensor_product_kernels.h:61
CUDAWrappers::FEEvaluation::distribute_local_to_global
void distribute_local_to_global(Number *dst) const
Definition: cuda_fe_evaluation.h:392
cuda_tensor_product_kernels.h
CUDAWrappers::FEEvaluation::JxW
Number * JxW
Definition: cuda_fe_evaluation.h:325
CUDAWrappers::FEEvaluation::n_q_points
static constexpr unsigned int n_q_points
Definition: cuda_fe_evaluation.h:123
CUDAWrappers::FEEvaluation::submit_gradient
void submit_gradient(const gradient_type &grad_in, const unsigned int q_point)
Definition: cuda_fe_evaluation.h:692
CUDAWrappers::FEEvaluation
Definition: cuda_fe_evaluation.h:92
CUDAWrappers::FEEvaluation::integrate
void integrate(const bool integrate_val, const bool integrate_grad)
Definition: cuda_fe_evaluation.h:457
tensor.h
CUDAWrappers::MatrixFree::Data
Definition: cuda_matrix_free.h:163
CUDAWrappers
Definition: cuda_size.h:23
Tensor
Definition: tensor.h:450
CUDAWrappers::FEEvaluation::gradients
Number * gradients[dim]
Definition: cuda_fe_evaluation.h:329
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
DEAL_II_DEPRECATED
#define DEAL_II_DEPRECATED
Definition: config.h:98
CUDAWrappers::internal::compute_index
unsigned int compute_index()
Definition: cuda_fe_evaluation.h:49
CUDAWrappers::FEEvaluation::submit_value
void submit_value(const value_type &val_in, const unsigned int q_point)
Definition: cuda_fe_evaluation.h:568
CUDAWrappers::FEEvaluation::get_value
value_type get_value() const
Definition: cuda_fe_evaluation.h:516
CUDAWrappers::FEEvaluation::get_gradient
gradient_type get_gradient() const
Definition: cuda_fe_evaluation.h:662
unsigned int
CUDAWrappers::FEEvaluation::n_cells
unsigned int n_cells
Definition: cuda_fe_evaluation.h:317
CUDAWrappers::FEEvaluation::use_coloring
const bool use_coloring
Definition: cuda_fe_evaluation.h:322
config.h
CUDAWrappers::FEEvaluation::tensor_dofs_per_cell
static constexpr unsigned int tensor_dofs_per_cell
Definition: cuda_fe_evaluation.h:129
internal
Definition: aligned_vector.h:369
cuda_hanging_nodes_internal.h
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12618
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
CUDAWrappers::FEEvaluation::apply_for_each_quad_point
void apply_for_each_quad_point(const Functor &func)
Definition: cuda_fe_evaluation.h:757
value_type
CUDAWrappers::FEEvaluation::constraint_mask
const unsigned int constraint_mask
Definition: cuda_fe_evaluation.h:320
CUDAWrappers::FEEvaluation::evaluate
void evaluate(const bool evaluate_val, const bool evaluate_grad)
Definition: cuda_fe_evaluation.h:418
cuda_atomic.h
CUDAWrappers::FEEvaluation::values
Number * values
Definition: cuda_fe_evaluation.h:328
CUDAWrappers::FEEvaluation::inv_jac
Number * inv_jac
Definition: cuda_fe_evaluation.h:324
Generated by   doxygen 1.8.17