Reference documentation for deal.II version 9.5.1
\(\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\}}\)
Loading...
Searching...
No Matches
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
23#include <deal.II/matrix_free/hanging_nodes_internal.h>
24
25#include <Kokkos_Macros.hpp>
26
27DEAL_II_NAMESPACE_OPEN
28namespace 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
55 is_constrained_dof_2d(
56 const ::internal::MatrixFreeFunctions::ConstraintKinds
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) !=
64 ::internal::MatrixFreeFunctions::ConstraintKinds::
65 unconstrained) ?
66 (y_idx == 0) :
67 (y_idx == fe_degree))) ||
68 ((direction == 1) &&
69 (((constraint_mask & ::internal::MatrixFreeFunctions::
70 ConstraintKinds::subcell_x) !=
71 ::internal::MatrixFreeFunctions::ConstraintKinds::
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
79 is_constrained_dof_3d(
80 const ::internal::MatrixFreeFunctions::ConstraintKinds
81 & constraint_mask,
82 const unsigned int x_idx,
83 const unsigned int y_idx,
84 const unsigned int z_idx,
85 const ::internal::MatrixFreeFunctions::ConstraintKinds face1_type,
86 const ::internal::MatrixFreeFunctions::ConstraintKinds face2_type,
87 const ::internal::MatrixFreeFunctions::ConstraintKinds face1,
88 const ::internal::MatrixFreeFunctions::ConstraintKinds face2,
89 const ::internal::MatrixFreeFunctions::ConstraintKinds edge)
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) !=
99 ::internal::MatrixFreeFunctions::
100 ConstraintKinds::unconstrained) ?
101 (face1_idx == 0) :
102 (face1_idx == fe_degree);
103 const bool on_face2 = ((constraint_mask & face2_type) !=
104 ::internal::MatrixFreeFunctions::
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
127 interpolate_boundary_2d(
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,
133 const ::internal::MatrixFreeFunctions::ConstraintKinds
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) ?
148 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
149 ::internal::MatrixFreeFunctions::ConstraintKinds::
150 unconstrained) |
151 ((direction == 1) ?
152 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
153 ::internal::MatrixFreeFunctions::ConstraintKinds::
154 unconstrained))) !=
155 ::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained;
156
157 Number tmp[n_q_points];
158 Kokkos::parallel_for(
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) ?
166 ::internal::MatrixFreeFunctions::ConstraintKinds::
167 subcell_x :
168 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y;
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) !=
184 ::internal::MatrixFreeFunctions::
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
248 interpolate_boundary_3d(
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,
254 const ::internal::MatrixFreeFunctions::ConstraintKinds
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) ?
266 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x :
267 (direction == 1) ?
268 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
269 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z;
270 const auto face1_type =
271 (direction == 0) ?
272 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
273 (direction == 1) ?
274 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
275 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x;
276 const auto face2_type =
277 (direction == 0) ?
278 ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
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) ?
287 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y :
288 (direction == 1) ?
289 ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
290 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x;
291 const auto face2 =
292 (direction == 0) ?
293 ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
294 (direction == 1) ?
295 ::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
296 ::internal::MatrixFreeFunctions::ConstraintKinds::face_y;
297 const auto edge =
298 (direction == 0) ?
299 ::internal::MatrixFreeFunctions::ConstraintKinds::edge_x :
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];
306 Kokkos::parallel_for(
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) !=
332 ::internal::MatrixFreeFunctions::
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
380 Kokkos::parallel_for(
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>
415 DEAL_II_HOST_DEVICE void
416 resolve_hanging_nodes(
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,
422 const ::internal::MatrixFreeFunctions::ConstraintKinds
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
463DEAL_II_NAMESPACE_CLOSE
464#endif
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
cuda_size.h
transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition derivative_form.h:586
hanging_nodes_internal.h
CUDAWrappers::internal::interpolate_boundary_3d
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)
Definition cuda_hanging_nodes_internal.h:248
CUDAWrappers::internal::index3
unsigned int index3(unsigned int i, unsigned int j, unsigned int k)
Definition cuda_hanging_nodes_internal.h:46
CUDAWrappers::internal::interpolate_boundary_2d
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)
Definition cuda_hanging_nodes_internal.h:127
CUDAWrappers::internal::is_constrained_dof_2d
bool is_constrained_dof_2d(const ::internal::MatrixFreeFunctions::ConstraintKinds &constraint_mask, const unsigned int x_idx, const unsigned int y_idx)
Definition cuda_hanging_nodes_internal.h:55
CUDAWrappers::internal::is_constrained_dof_3d
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)
Definition cuda_hanging_nodes_internal.h:79
CUDAWrappers::internal::index2
unsigned int index2(unsigned int i, unsigned int j)
Definition cuda_hanging_nodes_internal.h:37
CUDAWrappers::internal::resolve_hanging_nodes
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)
Definition cuda_hanging_nodes_internal.h:416
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition utilities.h:447
DEAL_II_HOST_DEVICE
#define DEAL_II_HOST_DEVICE
Definition numbers.h:35
Generated by doxygen 1.9.8