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