Reference documentation for deal.II version GIT 3c37cfefb2 2023-03-24 04:00:03+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 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 #ifdef DEAL_II_WITH_CUDA
22 
23 # include <deal.II/base/cuda_size.h>
24 
26 
28 namespace CUDAWrappers
29 {
30  namespace internal
31  {
32  __constant__ double
35 
36  //------------------------------------------------------------------------//
37  // Functions for resolving the hanging node constraints on the GPU //
38  //------------------------------------------------------------------------//
39  template <unsigned int size>
40  __device__ inline unsigned int
41  index2(unsigned int i, unsigned int j)
42  {
43  return i + size * j;
44  }
45 
46 
47 
48  template <unsigned int size>
49  __device__ inline unsigned int
50  index3(unsigned int i, unsigned int j, unsigned int k)
51  {
52  return i + size * j + size * size * k;
53  }
54 
55 
56 
57  template <unsigned int fe_degree,
58  unsigned int direction,
59  bool transpose,
60  typename Number>
61  __device__ inline void
64  constraint_mask,
65  Number *values)
66  {
67  const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
68  const unsigned int y_idx = threadIdx.y;
69 
70  const auto this_type =
71  (direction == 0) ?
74 
75  const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
76 
77  Number t = 0;
78  // Flag is true if dof is constrained for the given direction and the
79  // given face.
80  const bool constrained_face =
81  (constraint_mask &
82  (((direction == 0) ?
85  unconstrained) |
86  ((direction == 1) ?
89  unconstrained))) !=
90  ::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained;
91 
92  // Flag is true if for the given direction, the dof is constrained with
93  // the right type and is on the correct side (left (= 0) or right (=
94  // fe_degree))
95  const bool constrained_dof =
96  ((direction == 0) &&
97  (((constraint_mask & ::internal::MatrixFreeFunctions::
98  ConstraintKinds::subcell_y) !=
99  ::internal::MatrixFreeFunctions::ConstraintKinds::
100  unconstrained) ?
101  (y_idx == 0) :
102  (y_idx == fe_degree))) ||
103  ((direction == 1) &&
104  (((constraint_mask & ::internal::MatrixFreeFunctions::
105  ConstraintKinds::subcell_x) !=
107  unconstrained) ?
108  (x_idx == 0) :
109  (x_idx == fe_degree)));
110 
111  if (constrained_face && constrained_dof)
112  {
113  const bool type = (constraint_mask & this_type) !=
115  ConstraintKinds::unconstrained;
116 
117  if (type)
118  {
119  for (unsigned int i = 0; i <= fe_degree; ++i)
120  {
121  const unsigned int real_idx =
122  (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
123  index2<fe_degree + 1>(x_idx, i);
124 
125  const Number w =
126  transpose ?
127  constraint_weights[i * (fe_degree + 1) + interp_idx] :
128  constraint_weights[interp_idx * (fe_degree + 1) + i];
129  t += w * values[real_idx];
130  }
131  }
132  else
133  {
134  for (unsigned int i = 0; i <= fe_degree; ++i)
135  {
136  const unsigned int real_idx =
137  (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
138  index2<fe_degree + 1>(x_idx, i);
139 
140  const Number w =
141  transpose ?
142  constraint_weights[(fe_degree - i) * (fe_degree + 1) +
143  fe_degree - interp_idx] :
144  constraint_weights[(fe_degree - interp_idx) *
145  (fe_degree + 1) +
146  fe_degree - i];
147  t += w * values[real_idx];
148  }
149  }
150  }
151 
152  // The synchronization is done for all the threads in one block with
153  // each block being assigned to one element.
154  __syncthreads();
155  if (constrained_face && constrained_dof)
156  values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
157 
158  __syncthreads();
159  }
160 
161 
162 
163  template <unsigned int fe_degree,
164  unsigned int direction,
165  bool transpose,
166  typename Number>
167  __device__ inline void
170  constraint_mask,
171  Number *values)
172  {
173  const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
174  const unsigned int y_idx = threadIdx.y;
175  const unsigned int z_idx = threadIdx.z;
176 
177  const auto this_type =
178  (direction == 0) ?
180  (direction == 1) ?
181  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y :
182  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z;
183  const auto face1_type =
184  (direction == 0) ?
186  (direction == 1) ?
187  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_z :
188  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x;
189  const auto face2_type =
190  (direction == 0) ?
192  (direction == 1) ?
193  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_x :
194  ::internal::MatrixFreeFunctions::ConstraintKinds::subcell_y;
195 
196  // If computing in x-direction, need to match against face_y or
197  // face_z
198  const auto face1 =
199  (direction == 0) ?
201  (direction == 1) ?
202  ::internal::MatrixFreeFunctions::ConstraintKinds::face_z :
203  ::internal::MatrixFreeFunctions::ConstraintKinds::face_x;
204  const auto face2 =
205  (direction == 0) ?
207  (direction == 1) ?
208  ::internal::MatrixFreeFunctions::ConstraintKinds::face_x :
209  ::internal::MatrixFreeFunctions::ConstraintKinds::face_y;
210  const auto edge =
211  (direction == 0) ?
213  (direction == 1) ?
214  ::internal::MatrixFreeFunctions::ConstraintKinds::edge_y :
215  ::internal::MatrixFreeFunctions::ConstraintKinds::edge_z;
216  const auto constrained_face = constraint_mask & (face1 | face2 | edge);
217 
218  const unsigned int interp_idx = (direction == 0) ? x_idx :
219  (direction == 1) ? y_idx :
220  z_idx;
221  const unsigned int face1_idx = (direction == 0) ? y_idx :
222  (direction == 1) ? z_idx :
223  x_idx;
224  const unsigned int face2_idx = (direction == 0) ? z_idx :
225  (direction == 1) ? x_idx :
226  y_idx;
227 
228  Number t = 0;
229  const bool on_face1 = ((constraint_mask & face1_type) !=
231  ConstraintKinds::unconstrained) ?
232  (face1_idx == 0) :
233  (face1_idx == fe_degree);
234  const bool on_face2 = ((constraint_mask & face2_type) !=
236  ConstraintKinds::unconstrained) ?
237  (face2_idx == 0) :
238  (face2_idx == fe_degree);
239  const bool constrained_dof =
240  ((((constraint_mask & face1) != ::internal::MatrixFreeFunctions::
241  ConstraintKinds::unconstrained) &&
242  on_face1) ||
243  (((constraint_mask & face2) != ::internal::MatrixFreeFunctions::
244  ConstraintKinds::unconstrained) &&
245  on_face2) ||
246  (((constraint_mask & edge) != ::internal::MatrixFreeFunctions::
247  ConstraintKinds::unconstrained) &&
248  on_face1 && on_face2));
249 
250  if ((constrained_face != ::internal::MatrixFreeFunctions::
251  ConstraintKinds::unconstrained) &&
252  constrained_dof)
253  {
254  const bool type = (constraint_mask & this_type) !=
256  ConstraintKinds::unconstrained;
257  if (type)
258  {
259  for (unsigned int i = 0; i <= fe_degree; ++i)
260  {
261  const unsigned int real_idx =
262  (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
263  (direction == 1) ? index3<fe_degree + 1>(x_idx, i, z_idx) :
264  index3<fe_degree + 1>(x_idx, y_idx, i);
265 
266  const Number w =
267  transpose ?
268  constraint_weights[i * (fe_degree + 1) + interp_idx] :
269  constraint_weights[interp_idx * (fe_degree + 1) + i];
270  t += w * values[real_idx];
271  }
272  }
273  else
274  {
275  for (unsigned int i = 0; i <= fe_degree; ++i)
276  {
277  const unsigned int real_idx =
278  (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
279  (direction == 1) ? index3<fe_degree + 1>(x_idx, i, z_idx) :
280  index3<fe_degree + 1>(x_idx, y_idx, i);
281 
282  const Number w =
283  transpose ?
284  constraint_weights[(fe_degree - i) * (fe_degree + 1) +
285  fe_degree - interp_idx] :
286  constraint_weights[(fe_degree - interp_idx) *
287  (fe_degree + 1) +
288  fe_degree - i];
289  t += w * values[real_idx];
290  }
291  }
292  }
293 
294  // The synchronization is done for all the threads in one block with
295  // each block being assigned to one element.
296  __syncthreads();
297 
298  if ((constrained_face != ::internal::MatrixFreeFunctions::
299  ConstraintKinds::unconstrained) &&
300  constrained_dof)
301  values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
302 
303  __syncthreads();
304  }
305 
306 
307 
315  template <int dim, int fe_degree, bool transpose, typename Number>
316  __device__ void
319  constraint_mask,
320  Number *values)
321  {
322  if (dim == 2)
323  {
324  interpolate_boundary_2d<fe_degree, 0, transpose>(constraint_mask,
325  values);
326 
327  interpolate_boundary_2d<fe_degree, 1, transpose>(constraint_mask,
328  values);
329  }
330  else if (dim == 3)
331  {
332  // Interpolate y and z faces (x-direction)
333  interpolate_boundary_3d<fe_degree, 0, transpose>(constraint_mask,
334  values);
335  // Interpolate x and z faces (y-direction)
336  interpolate_boundary_3d<fe_degree, 1, transpose>(constraint_mask,
337  values);
338  // Interpolate x and y faces (z-direction)
339  interpolate_boundary_3d<fe_degree, 2, transpose>(constraint_mask,
340  values);
341  }
342  }
343  } // namespace internal
344 } // namespace CUDAWrappers
345 
347 #endif
348 
349 #endif
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
void interpolate_boundary_2d(const ::internal::MatrixFreeFunctions::ConstraintKinds constraint_mask, Number *values)
void resolve_hanging_nodes(const ::internal::MatrixFreeFunctions::ConstraintKinds constraint_mask, Number *values)
void interpolate_boundary_3d(const ::internal::MatrixFreeFunctions::ConstraintKinds constraint_mask, Number *values)
unsigned int index3(unsigned int i, unsigned int j, unsigned int k)
__constant__ double constraint_weights[(CUDAWrappers::mf_max_elem_degree+1) *(CUDAWrappers::mf_max_elem_degree+1)]
unsigned int index2(unsigned int i, unsigned int j)
constexpr unsigned int mf_max_elem_degree
Definition: cuda_size.h:47
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)