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