Reference documentation for deal.II version GIT c00406db16 2023-09-28 18:00: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_hanging_nodes_internal.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2021 - 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_cuda_hanging_nodes_internal_h
17 #define dealii_cuda_hanging_nodes_internal_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/cuda_size.h>
22 
24 
25 #include <Kokkos_Macros.hpp>
26 
28 namespace CUDAWrappers
29 {
30  namespace internal
31  {
32  //------------------------------------------------------------------------//
33  // Functions for resolving the hanging node constraints on the GPU //
34  //------------------------------------------------------------------------//
35  template <unsigned int size>
36  DEAL_II_HOST_DEVICE inline unsigned int
37  index2(unsigned int i, unsigned int j)
38  {
39  return i + size * j;
40  }
41 
42 
43 
44  template <unsigned int size>
45  DEAL_II_HOST_DEVICE inline unsigned int
46  index3(unsigned int i, unsigned int j, unsigned int k)
47  {
48  return i + size * j + size * size * k;
49  }
50 
51 
52 
53  template <unsigned int fe_degree, unsigned int direction>
54  DEAL_II_HOST_DEVICE inline bool
57  &constraint_mask,
58  const unsigned int x_idx,
59  const unsigned int y_idx)
60  {
61  return ((direction == 0) &&
62  (((constraint_mask & ::internal::MatrixFreeFunctions::
63  ConstraintKinds::subcell_y) !=
65  unconstrained) ?
66  (y_idx == 0) :
67  (y_idx == fe_degree))) ||
68  ((direction == 1) &&
69  (((constraint_mask & ::internal::MatrixFreeFunctions::
70  ConstraintKinds::subcell_x) !=
72  unconstrained) ?
73  (x_idx == 0) :
74  (x_idx == fe_degree)));
75  }
76 
77  template <unsigned int fe_degree, unsigned int direction>
78  DEAL_II_HOST_DEVICE inline bool
81  &constraint_mask,
82  const unsigned int x_idx,
83  const unsigned int y_idx,
84  const unsigned int z_idx,
90  {
91  const unsigned int face1_idx = (direction == 0) ? y_idx :
92  (direction == 1) ? z_idx :
93  x_idx;
94  const unsigned int face2_idx = (direction == 0) ? z_idx :
95  (direction == 1) ? x_idx :
96  y_idx;
97 
98  const bool on_face1 = ((constraint_mask & face1_type) !=
100  ConstraintKinds::unconstrained) ?
101  (face1_idx == 0) :
102  (face1_idx == fe_degree);
103  const bool on_face2 = ((constraint_mask & face2_type) !=
105  ConstraintKinds::unconstrained) ?
106  (face2_idx == 0) :
107  (face2_idx == fe_degree);
108  return (
109  (((constraint_mask & face1) != ::internal::MatrixFreeFunctions::
110  ConstraintKinds::unconstrained) &&
111  on_face1) ||
112  (((constraint_mask & face2) != ::internal::MatrixFreeFunctions::
113  ConstraintKinds::unconstrained) &&
114  on_face2) ||
115  (((constraint_mask & edge) != ::internal::MatrixFreeFunctions::
116  ConstraintKinds::unconstrained) &&
117  on_face1 && on_face2));
118  }
119 
120 
121 
122  template <unsigned int fe_degree,
123  unsigned int direction,
124  bool transpose,
125  typename Number>
126  DEAL_II_HOST_DEVICE inline void
128  const Kokkos::TeamPolicy<
129  MemorySpace::Default::kokkos_space::execution_space>::member_type
130  &team_member,
131  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
132  constraint_weights,
134  &constraint_mask,
135  Kokkos::View<Number *,
136  MemorySpace::Default::kokkos_space::execution_space::
137  scratch_memory_space,
138  Kokkos::MemoryTraits<Kokkos::Unmanaged>> values)
139  {
140  constexpr unsigned int n_q_points_1d = fe_degree + 1;
141  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 2);
142 
143  // Flag is true if dof is constrained for the given direction and the
144  // given face.
145  const bool constrained_face =
146  (constraint_mask &
147  (((direction == 0) ?
150  unconstrained) |
151  ((direction == 1) ?
154  unconstrained))) !=
155  ::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained;
156 
157  Number tmp[n_q_points];
159  Kokkos::TeamThreadRange(team_member, n_q_points),
160  [&](const int &q_point) {
161  const unsigned int x_idx = q_point % n_q_points_1d;
162  const unsigned int y_idx = q_point / n_q_points_1d;
163 
164  const auto this_type =
165  (direction == 0) ?
167  subcell_x :
169 
170  const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
171  tmp[q_point] = 0;
172 
173  // Flag is true if for the given direction, the dof is constrained
174  // with the right type and is on the correct side (left (= 0) or right
175  // (= fe_degree))
176  const bool constrained_dof =
177  is_constrained_dof_2d<fe_degree, direction>(constraint_mask,
178  x_idx,
179  y_idx);
180 
181  if (constrained_face && constrained_dof)
182  {
183  const bool type = (constraint_mask & this_type) !=
185  ConstraintKinds::unconstrained;
186 
187  if (type)
188  {
189  for (unsigned int i = 0; i <= fe_degree; ++i)
190  {
191  const unsigned int real_idx =
192  (direction == 0) ? index2<n_q_points_1d>(i, y_idx) :
193  index2<n_q_points_1d>(x_idx, i);
194 
195  const Number w =
196  transpose ?
197  constraint_weights[i * n_q_points_1d + interp_idx] :
198  constraint_weights[interp_idx * n_q_points_1d + i];
199  tmp[q_point] += w * values[real_idx];
200  }
201  }
202  else
203  {
204  for (unsigned int i = 0; i <= fe_degree; ++i)
205  {
206  const unsigned int real_idx =
207  (direction == 0) ? index2<n_q_points_1d>(i, y_idx) :
208  index2<n_q_points_1d>(x_idx, i);
209 
210  const Number w =
211  transpose ?
212  constraint_weights[(fe_degree - i) * n_q_points_1d +
213  fe_degree - interp_idx] :
214  constraint_weights[(fe_degree - interp_idx) *
215  n_q_points_1d +
216  fe_degree - i];
217  tmp[q_point] += w * values[real_idx];
218  }
219  }
220  }
221  });
222 
223  // The synchronization is done for all the threads in one team with
224  // each team being assigned to one element.
225  team_member.team_barrier();
226  Kokkos::parallel_for(Kokkos::TeamThreadRange(team_member, n_q_points),
227  [&](const int &q_point) {
228  const unsigned int x_idx = q_point % n_q_points_1d;
229  const unsigned int y_idx = q_point / n_q_points_1d;
230  const bool constrained_dof =
231  is_constrained_dof_2d<fe_degree, direction>(
232  constraint_mask, x_idx, y_idx);
233  if (constrained_face && constrained_dof)
234  values[index2<fe_degree + 1>(x_idx, y_idx)] =
235  tmp[q_point];
236  });
237 
238  team_member.team_barrier();
239  }
240 
241 
242 
243  template <unsigned int fe_degree,
244  unsigned int direction,
245  bool transpose,
246  typename Number>
247  DEAL_II_HOST_DEVICE inline void
249  const Kokkos::TeamPolicy<
250  MemorySpace::Default::kokkos_space::execution_space>::member_type
251  &team_member,
252  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
253  constraint_weights,
255  constraint_mask,
256  Kokkos::View<Number *,
257  MemorySpace::Default::kokkos_space::execution_space::
258  scratch_memory_space,
259  Kokkos::MemoryTraits<Kokkos::Unmanaged>> values)
260  {
261  constexpr unsigned int n_q_points_1d = fe_degree + 1;
262  constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, 3);
263 
264  const auto this_type =
265  (direction == 0) ?
267  (direction == 1) ?
268  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
269  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z;
270  const auto face1_type =
271  (direction == 0) ?
273  (direction == 1) ?
274  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
275  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x;
276  const auto face2_type =
277  (direction == 0) ?
279  (direction == 1) ?
280  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x :
281  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y;
282 
283  // If computing in x-direction, need to match against face_y or
284  // face_z
285  const auto face1 =
286  (direction == 0) ?
288  (direction == 1) ?
289  ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
290  ::internal::MatrixFreeFunctions::ConstraintKinds::face_x;
291  const auto face2 =
292  (direction == 0) ?
294  (direction == 1) ?
295  ::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
296  ::internal::MatrixFreeFunctions::ConstraintKinds::face_y;
297  const auto edge =
298  (direction == 0) ?
300  (direction == 1) ?
301  ::internal::MatrixFreeFunctions::ConstraintKinds::edge_y :
302  ::internal::MatrixFreeFunctions::ConstraintKinds::edge_z;
303  const auto constrained_face = constraint_mask & (face1 | face2 | edge);
304 
305  Number tmp[n_q_points];
307  Kokkos::TeamThreadRange(team_member, n_q_points),
308  [&](const int &q_point) {
309  const unsigned int x_idx = q_point % n_q_points_1d;
310  const unsigned int y_idx = (q_point / n_q_points_1d) % n_q_points_1d;
311  const unsigned int z_idx = q_point / (n_q_points_1d * n_q_points_1d);
312 
313  const unsigned int interp_idx = (direction == 0) ? x_idx :
314  (direction == 1) ? y_idx :
315  z_idx;
316  const bool constrained_dof =
317  is_constrained_dof_3d<fe_degree, direction>(constraint_mask,
318  x_idx,
319  y_idx,
320  z_idx,
321  face1_type,
322  face2_type,
323  face1,
324  face2,
325  edge);
326  tmp[q_point] = 0;
327  if ((constrained_face != ::internal::MatrixFreeFunctions::
328  ConstraintKinds::unconstrained) &&
329  constrained_dof)
330  {
331  const bool type = (constraint_mask & this_type) !=
333  ConstraintKinds::unconstrained;
334  if (type)
335  {
336  for (unsigned int i = 0; i <= fe_degree; ++i)
337  {
338  const unsigned int real_idx =
339  (direction == 0) ?
340  index3<fe_degree + 1>(i, y_idx, z_idx) :
341  (direction == 1) ?
342  index3<fe_degree + 1>(x_idx, i, z_idx) :
343  index3<fe_degree + 1>(x_idx, y_idx, i);
344 
345  const Number w =
346  transpose ?
347  constraint_weights[i * n_q_points_1d + interp_idx] :
348  constraint_weights[interp_idx * n_q_points_1d + i];
349  tmp[q_point] += w * values[real_idx];
350  }
351  }
352  else
353  {
354  for (unsigned int i = 0; i <= fe_degree; ++i)
355  {
356  const unsigned int real_idx =
357  (direction == 0) ?
358  index3<n_q_points_1d>(i, y_idx, z_idx) :
359  (direction == 1) ?
360  index3<n_q_points_1d>(x_idx, i, z_idx) :
361  index3<n_q_points_1d>(x_idx, y_idx, i);
362 
363  const Number w =
364  transpose ?
365  constraint_weights[(fe_degree - i) * n_q_points_1d +
366  fe_degree - interp_idx] :
367  constraint_weights[(fe_degree - interp_idx) *
368  n_q_points_1d +
369  fe_degree - i];
370  tmp[q_point] += w * values[real_idx];
371  }
372  }
373  }
374  });
375 
376  // The synchronization is done for all the threads in one team with
377  // each team being assigned to one element.
378  team_member.team_barrier();
379 
381  Kokkos::TeamThreadRange(team_member, n_q_points),
382  [&](const int &q_point) {
383  const unsigned int x_idx = q_point % n_q_points_1d;
384  const unsigned int y_idx = (q_point / n_q_points_1d) % n_q_points_1d;
385  const unsigned int z_idx = q_point / (n_q_points_1d * n_q_points_1d);
386  const bool constrained_dof =
387  is_constrained_dof_3d<fe_degree, direction>(constraint_mask,
388  x_idx,
389  y_idx,
390  z_idx,
391  face1_type,
392  face2_type,
393  face1,
394  face2,
395  edge);
396  if ((constrained_face != ::internal::MatrixFreeFunctions::
397  ConstraintKinds::unconstrained) &&
398  constrained_dof)
399  values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = tmp[q_point];
400  });
401 
402  team_member.team_barrier();
403  }
404 
405 
406 
414  template <int dim, int fe_degree, bool transpose, typename Number>
417  const Kokkos::TeamPolicy<
418  MemorySpace::Default::kokkos_space::execution_space>::member_type
419  &team_member,
420  Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
421  constraint_weights,
423  constraint_mask,
424  Kokkos::View<Number *,
425  MemorySpace::Default::kokkos_space::execution_space::
426  scratch_memory_space,
427  Kokkos::MemoryTraits<Kokkos::Unmanaged>> values)
428  {
429  if (dim == 2)
430  {
431  interpolate_boundary_2d<fe_degree, 0, transpose>(team_member,
432  constraint_weights,
433  constraint_mask,
434  values);
435 
436  interpolate_boundary_2d<fe_degree, 1, transpose>(team_member,
437  constraint_weights,
438  constraint_mask,
439  values);
440  }
441  else if (dim == 3)
442  {
443  // Interpolate y and z faces (x-direction)
444  interpolate_boundary_3d<fe_degree, 0, transpose>(team_member,
445  constraint_weights,
446  constraint_mask,
447  values);
448  // Interpolate x and z faces (y-direction)
449  interpolate_boundary_3d<fe_degree, 1, transpose>(team_member,
450  constraint_weights,
451  constraint_mask,
452  values);
453  // Interpolate x and y faces (z-direction)
454  interpolate_boundary_3d<fe_degree, 2, transpose>(team_member,
455  constraint_weights,
456  constraint_mask,
457  values);
458  }
459  }
460  } // namespace internal
461 } // namespace CUDAWrappers
462 
464 #endif
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
void interpolate_boundary_3d(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights, const ::internal::MatrixFreeFunctions::ConstraintKinds constraint_mask, Kokkos::View< Number *, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged >> values)
unsigned int index3(unsigned int i, unsigned int j, unsigned int k)
void interpolate_boundary_2d(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights, const ::internal::MatrixFreeFunctions::ConstraintKinds &constraint_mask, Kokkos::View< Number *, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged >> values)
void resolve_hanging_nodes(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, Kokkos::View< Number *, MemorySpace::Default::kokkos_space > constraint_weights, const ::internal::MatrixFreeFunctions::ConstraintKinds constraint_mask, Kokkos::View< Number *, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged >> values)
bool is_constrained_dof_2d(const ::internal::MatrixFreeFunctions::ConstraintKinds &constraint_mask, const unsigned int x_idx, const unsigned int y_idx)
bool is_constrained_dof_3d(const ::internal::MatrixFreeFunctions::ConstraintKinds &constraint_mask, const unsigned int x_idx, const unsigned int y_idx, const unsigned int z_idx, const ::internal::MatrixFreeFunctions::ConstraintKinds face1_type, const ::internal::MatrixFreeFunctions::ConstraintKinds face2_type, const ::internal::MatrixFreeFunctions::ConstraintKinds face1, const ::internal::MatrixFreeFunctions::ConstraintKinds face2, const ::internal::MatrixFreeFunctions::ConstraintKinds edge)
unsigned int index2(unsigned int i, unsigned int j)
typename ::internal::FEInterfaceViews::ViewType< dim, spacedim, Extractor >::type View
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
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