Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10: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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
portable_evaluation_kernels.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2017 - 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
16#ifndef dealii__evaluation_kernels_h
17#define dealii__evaluation_kernels_h
18
19#include <deal.II/base/config.h>
20
22
26
27#include <Kokkos_Core.hpp>
28
29
31
32
33namespace Portable
34{
35 namespace internal
36 {
44 constexpr bool
45 use_collocation_evaluation(const unsigned int fe_degree,
46 const unsigned int n_q_points_1d)
47 {
48 // TODO: are the conditions suit for GPU parallelization?
49 return (n_q_points_1d > fe_degree) && (n_q_points_1d < 200) &&
50 (n_q_points_1d <= 3 * fe_degree / 2 + 1);
51 }
52
53
54
64 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
66 {
67 using TeamHandle = Kokkos::TeamPolicy<
68 MemorySpace::Default::kokkos_space::execution_space>::member_type;
69 using SharedView = Kokkos::View<Number *,
70 MemorySpace::Default::kokkos_space::
71 execution_space::scratch_memory_space,
72 Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
73
74 DEAL_II_HOST_DEVICE static void
75 evaluate(const unsigned int n_components,
76 const EvaluationFlags::EvaluationFlags evaluation_flag,
77 const typename MatrixFree<dim, Number>::Data *data)
78 {
79 if (evaluation_flag == EvaluationFlags::nothing)
80 return;
81
82 // the evaluator does not need temporary storage since no in-place
83 // operation takes place in this function
84 auto scratch_for_eval = Kokkos::subview(data->shared_data->scratch_pad,
85 Kokkos::make_pair(0, 0));
87 dim,
88 fe_degree + 1,
89 n_q_points_1d,
90 Number>
91 eval(data->team_member,
92 data->precomputed_data->shape_values,
93 data->precomputed_data->shape_gradients,
94 data->precomputed_data->co_shape_gradients,
95 scratch_for_eval);
96
97 for (unsigned int c = 0; c < n_components; ++c)
98 {
99 auto u = Kokkos::subview(data->shared_data->values, Kokkos::ALL, c);
100 auto grad_u = Kokkos::subview(data->shared_data->gradients,
101 Kokkos::ALL,
102 Kokkos::ALL,
103 c);
104
105 if constexpr (dim == 1)
106 {
107 auto temp =
108 Kokkos::subview(data->shared_data->scratch_pad,
109 Kokkos::make_pair(0, n_q_points_1d));
110
111 if (evaluation_flag & EvaluationFlags::gradients)
112 eval.template gradients<0, true, false, false>(
113 u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
114 if (evaluation_flag & EvaluationFlags::values)
115 {
116 eval.template values<0, true, false, false>(u, temp);
117 populate_view<false>(data->team_member,
118 u,
119 temp,
120 n_q_points_1d);
121 }
122 }
123 else if constexpr (dim == 2)
124 {
125 constexpr int temp_size = (fe_degree + 1) * n_q_points_1d;
126 auto temp = Kokkos::subview(data->shared_data->scratch_pad,
127 Kokkos::make_pair(0, temp_size));
128
129 // grad x
130 if (evaluation_flag & EvaluationFlags::gradients)
131 {
132 eval.template gradients<0, true, false, false>(u, temp);
133 eval.template values<1, true, false, false>(
134 temp, Kokkos::subview(grad_u, Kokkos::ALL, 0));
135 }
136
137 // grad y
138 eval.template values<0, true, false, false>(u, temp);
139 if (evaluation_flag & EvaluationFlags::gradients)
140 eval.template gradients<1, true, false, false>(
141 temp, Kokkos::subview(grad_u, Kokkos::ALL, 1));
142
143 // val: can use values applied in x
144 if (evaluation_flag & EvaluationFlags::values)
145 eval.template values<1, true, false, false>(temp, u);
146 }
147 else if constexpr (dim == 3)
148 {
149 constexpr int temp1_size = Utilities::pow(fe_degree + 1, 2) *
150 n_q_points_1d,
151 temp2_size = Utilities::pow(n_q_points_1d, 2) *
152 (fe_degree + 1);
153
154 auto temp1 = Kokkos::subview(data->shared_data->scratch_pad,
155 Kokkos::make_pair(0, temp1_size));
156 auto temp2 =
157 Kokkos::subview(data->shared_data->scratch_pad,
158 Kokkos::make_pair(temp1_size,
159 temp1_size + temp2_size));
160
161 if (evaluation_flag & EvaluationFlags::gradients)
162 {
163 // grad x
164 eval.template gradients<0, true, false, false>(u, temp1);
165 eval.template values<1, true, false, false>(temp1, temp2);
166 eval.template values<2, true, false, false>(
167 temp2, Kokkos::subview(grad_u, Kokkos::ALL, 0));
168 }
169
170 // grad y
171 eval.template values<0, true, false, false>(u, temp1);
172 if (evaluation_flag & EvaluationFlags::gradients)
173 {
174 eval.template gradients<1, true, false, false>(temp1,
175 temp2);
176 eval.template values<2, true, false, false>(
177 temp2, Kokkos::subview(grad_u, Kokkos::ALL, 1));
178 }
179
180 // grad z: can use the values applied in x direction stored
181 // in temp1
182 eval.template values<1, true, false, false>(temp1, temp2);
183 if (evaluation_flag & EvaluationFlags::gradients)
184 eval.template gradients<2, true, false, false>(
185 temp2, Kokkos::subview(grad_u, Kokkos::ALL, 2));
186
187 // val: can use the values applied in x & y direction stored
188 // in temp2
189 if (evaluation_flag & EvaluationFlags::values)
190 eval.template values<2, true, false, false>(temp2, u);
191 }
192 else
193 Assert(false, ExcMessage("dim must not exceed 3!"));
194 }
195 }
196
197
198
199 DEAL_II_HOST_DEVICE static void
200 integrate(const unsigned int n_components,
201 const EvaluationFlags::EvaluationFlags integration_flag,
202 const typename MatrixFree<dim, Number>::Data *data)
203 {
204 if (integration_flag == EvaluationFlags::nothing)
205 return;
206
207 // the evaluator does not need temporary storage since no in-place
208 // operation takes place in this function
209 auto scratch_for_eval = Kokkos::subview(data->shared_data->scratch_pad,
210 Kokkos::make_pair(0, 0));
212 dim,
213 fe_degree + 1,
214 n_q_points_1d,
215 Number>
216 eval(data->team_member,
217 data->precomputed_data->shape_values,
218 data->precomputed_data->shape_gradients,
219 data->precomputed_data->co_shape_gradients,
220 scratch_for_eval);
221
222 for (unsigned int c = 0; c < n_components; ++c)
223 {
224 auto u = Kokkos::subview(data->shared_data->values, Kokkos::ALL, c);
225 auto grad_u = Kokkos::subview(data->shared_data->gradients,
226 Kokkos::ALL,
227 Kokkos::ALL,
228 c);
229
230 if constexpr (dim == 1)
231 {
232 auto temp =
233 Kokkos::subview(data->shared_data->scratch_pad,
234 Kokkos::make_pair(0, fe_degree + 1));
235
236 if ((integration_flag & EvaluationFlags::values) &&
237 !(integration_flag & EvaluationFlags::gradients))
238 {
239 eval.template values<0, false, false, false>(u, temp);
240 populate_view<false>(data->team_member,
241 u,
242 temp,
243 fe_degree + 1);
244 }
245 if (integration_flag & EvaluationFlags::gradients)
246 {
247 if (integration_flag & EvaluationFlags::values)
248 {
249 eval.template values<0, false, false, false>(u, temp);
250 eval.template gradients<0, false, true, false>(
251 Kokkos::subview(grad_u, Kokkos::ALL, 0), temp);
252 populate_view<false>(data->team_member,
253 u,
254 temp,
255 fe_degree + 1);
256 }
257 else
258 eval.template gradients<0, false, false, false>(
259 Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
260 }
261 }
262 else if constexpr (dim == 2)
263 {
264 constexpr int temp_size = (fe_degree + 1) * n_q_points_1d;
265 auto temp = Kokkos::subview(data->shared_data->scratch_pad,
266 Kokkos::make_pair(0, temp_size));
267
268 if ((integration_flag & EvaluationFlags::values) &&
269 !(integration_flag & EvaluationFlags::gradients))
270 {
271 eval.template values<1, false, false, false>(u, temp);
272 eval.template values<0, false, false, false>(temp, u);
273 }
274 if (integration_flag & EvaluationFlags::gradients)
275 {
276 eval.template gradients<1, false, false, false>(
277 Kokkos::subview(grad_u, Kokkos::ALL, 1), temp);
278 if (integration_flag & EvaluationFlags::values)
279 eval.template values<1, false, true, false>(u, temp);
280 eval.template values<0, false, false, false>(temp, u);
281 eval.template values<1, false, false, false>(
282 Kokkos::subview(grad_u, Kokkos::ALL, 0), temp);
283 eval.template gradients<0, false, true, false>(temp, u);
284 }
285 }
286 else if constexpr (dim == 3)
287 {
288 constexpr int temp1_size = Utilities::pow(n_q_points_1d, 2) *
289 (fe_degree + 1),
290 temp2_size = Utilities::pow(fe_degree + 1, 2) *
291 n_q_points_1d;
292
293 auto temp1 = Kokkos::subview(data->shared_data->scratch_pad,
294 Kokkos::make_pair(0, temp1_size));
295 auto temp2 =
296 Kokkos::subview(data->shared_data->scratch_pad,
297 Kokkos::make_pair(temp1_size,
298 temp1_size + temp2_size));
299
300 if ((integration_flag & EvaluationFlags::values) &&
301 !(integration_flag & EvaluationFlags::gradients))
302 {
303 eval.template values<2, false, false, false>(u, temp1);
304 eval.template values<1, false, false, false>(temp1, temp2);
305 eval.template values<0, false, false, false>(temp2, u);
306 }
307 if (integration_flag & EvaluationFlags::gradients)
308 {
309 eval.template gradients<2, false, false, false>(
310 Kokkos::subview(grad_u, Kokkos::ALL, 2), temp1);
311 if (integration_flag & EvaluationFlags::values)
312 eval.template values<2, false, true, false>(u, temp1);
313 eval.template values<1, false, false, false>(temp1, temp2);
314 eval.template values<2, false, false, false>(
315 Kokkos::subview(grad_u, Kokkos::ALL, 1), temp1);
316 eval.template gradients<1, false, true, false>(temp1,
317 temp2);
318 eval.template values<0, false, false, false>(temp2, u);
319 eval.template values<2, false, false, false>(
320 Kokkos::subview(grad_u, Kokkos::ALL, 0), temp1);
321 eval.template values<1, false, false, false>(temp1, temp2);
322 eval.template gradients<0, false, true, false>(temp2, u);
323 }
324 }
325 else
326 Assert(false, ExcMessage("dim must not exceed 3!"));
327 }
328 }
329 };
330
331
332
345 template <int dim, int fe_degree, typename Number>
347 {
348 DEAL_II_HOST_DEVICE static void
349 evaluate(const unsigned int n_components,
350 const EvaluationFlags::EvaluationFlags evaluation_flag,
351 const typename MatrixFree<dim, Number>::Data *data)
352 {
353 // since the dof values have already been stored in
354 // shared_data->values, there is nothing to do if the gradients are
355 // not required
356 if (!(evaluation_flag & EvaluationFlags::gradients))
357 return;
358
359 constexpr int n_points = Utilities::pow(fe_degree + 1, dim);
360 auto scratch_for_eval = Kokkos::subview(data->shared_data->scratch_pad,
361 Kokkos::make_pair(0, n_points));
362
364 dim,
365 fe_degree + 1,
366 fe_degree + 1,
367 Number>
368 eval(data->team_member,
369 data->precomputed_data->shape_values,
370 data->precomputed_data->shape_gradients,
371 data->precomputed_data->co_shape_gradients,
372 scratch_for_eval);
373
374 for (unsigned int c = 0; c < n_components; ++c)
375 {
376 auto u = Kokkos::subview(data->shared_data->values, Kokkos::ALL, c);
377 auto grad_u = Kokkos::subview(data->shared_data->gradients,
378 Kokkos::ALL,
379 Kokkos::ALL,
380 c);
381
382 eval.template co_gradients<0, true, false, false>(
383 u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
384 if constexpr (dim > 1)
385 eval.template co_gradients<1, true, false, false>(
386 u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
387 if constexpr (dim > 2)
388 eval.template co_gradients<2, true, false, false>(
389 u, Kokkos::subview(grad_u, Kokkos::ALL, 2));
390 }
391 }
392
393
394 DEAL_II_HOST_DEVICE static void
395 integrate(const unsigned int n_components,
396 const EvaluationFlags::EvaluationFlags integration_flag,
397 const typename MatrixFree<dim, Number>::Data *data)
398 {
399 // since the quad values have already been stored in
400 // shared_data->values, there is nothing to do if the gradients are
401 // not required
402 if (!(integration_flag & EvaluationFlags::gradients))
403 return;
404
405 constexpr int n_points = Utilities::pow(fe_degree + 1, dim);
406 auto scratch_for_eval = Kokkos::subview(data->shared_data->scratch_pad,
407 Kokkos::make_pair(0, n_points));
408
410 dim,
411 fe_degree + 1,
412 fe_degree + 1,
413 Number>
414 eval(data->team_member,
415 data->precomputed_data->shape_values,
416 data->precomputed_data->shape_gradients,
417 data->precomputed_data->co_shape_gradients,
418 scratch_for_eval);
419
420 for (unsigned int c = 0; c < n_components; ++c)
421 {
422 auto u = Kokkos::subview(data->shared_data->values, Kokkos::ALL, c);
423 auto grad_u = Kokkos::subview(data->shared_data->gradients,
424 Kokkos::ALL,
425 Kokkos::ALL,
426 c);
427
428 if constexpr (dim == 1)
429 {
430 if (integration_flag & EvaluationFlags::values)
431 eval.template co_gradients<0, false, true, false>(
432 Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
433 else
434 eval.template co_gradients<2, false, false, false>(
435 Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
436 }
437 else if constexpr (dim == 2)
438 {
439 if (integration_flag & EvaluationFlags::values)
440 eval.template co_gradients<1, false, true, false>(
441 Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
442 else
443 eval.template co_gradients<1, false, false, false>(
444 Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
445 eval.template co_gradients<0, false, true, false>(
446 Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
447 }
448 else if constexpr (dim == 3)
449 {
450 if (integration_flag & EvaluationFlags::values)
451 eval.template co_gradients<2, false, true, false>(
452 Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
453 else
454 eval.template co_gradients<2, false, false, false>(
455 Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
456 eval.template co_gradients<1, false, true, false>(
457 Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
458 eval.template co_gradients<0, false, true, false>(
459 Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
460 }
461 else
462 Assert(false, ExcMessage("dim must not exceed 3!"));
463 }
464 }
465 };
466
467
468
479 template <int dim, int fe_degree, int n_q_points_1d, typename Number>
481 {
482 DEAL_II_HOST_DEVICE static void
483 evaluate(const unsigned int n_components,
484 const EvaluationFlags::EvaluationFlags evaluation_flag,
485 const typename MatrixFree<dim, Number>::Data *data)
486 {
487 constexpr int scratch_size = Utilities::pow(n_q_points_1d, dim);
488 auto scratch_for_eval =
489 Kokkos::subview(data->shared_data->scratch_pad,
490 Kokkos::make_pair(0, scratch_size));
491
493 dim,
494 fe_degree + 1,
495 n_q_points_1d,
496 Number>
497 eval(data->team_member,
498 data->precomputed_data->shape_values,
499 data->precomputed_data->shape_gradients,
500 data->precomputed_data->co_shape_gradients,
501 scratch_for_eval);
502
503 for (unsigned int c = 0; c < n_components; ++c)
504 {
505 auto u = Kokkos::subview(data->shared_data->values, Kokkos::ALL, c);
506 auto grad_u = Kokkos::subview(data->shared_data->gradients,
507 Kokkos::ALL,
508 Kokkos::ALL,
509 c);
510
511 eval.template values<0, true, false, true>(u, u);
512 if constexpr (dim > 1)
513 eval.template values<1, true, false, true>(u, u);
514 if constexpr (dim > 2)
515 eval.template values<2, true, false, true>(u, u);
516
517 if (evaluation_flag & EvaluationFlags::gradients)
518 {
519 eval.template co_gradients<0, true, false, false>(
520 u, Kokkos::subview(grad_u, Kokkos::ALL, 0));
521 if constexpr (dim > 1)
522 eval.template co_gradients<1, true, false, false>(
523 u, Kokkos::subview(grad_u, Kokkos::ALL, 1));
524 if constexpr (dim > 2)
525 eval.template co_gradients<2, true, false, false>(
526 u, Kokkos::subview(grad_u, Kokkos::ALL, 2));
527 }
528 }
529 }
530
531
532 DEAL_II_HOST_DEVICE static void
533 integrate(const unsigned int n_components,
534 const EvaluationFlags::EvaluationFlags integration_flag,
535 const typename MatrixFree<dim, Number>::Data *data)
536 {
537 constexpr int scratch_size = Utilities::pow(n_q_points_1d, dim);
538 auto scratch_for_eval =
539 Kokkos::subview(data->shared_data->scratch_pad,
540 Kokkos::make_pair(0, scratch_size));
541
543 dim,
544 fe_degree + 1,
545 n_q_points_1d,
546 Number>
547 eval(data->team_member,
548 data->precomputed_data->shape_values,
549 data->precomputed_data->shape_gradients,
550 data->precomputed_data->co_shape_gradients,
551 scratch_for_eval);
552
553 for (unsigned int c = 0; c < n_components; ++c)
554 {
555 auto u = Kokkos::subview(data->shared_data->values, Kokkos::ALL, c);
556 auto grad_u = Kokkos::subview(data->shared_data->gradients,
557 Kokkos::ALL,
558 Kokkos::ALL,
559 c);
560
561 // apply derivatives in collocation space
562 if (integration_flag & EvaluationFlags::gradients)
563 {
564 if constexpr (dim == 1)
565 {
566 if (integration_flag & EvaluationFlags::values)
567 eval.template co_gradients<0, false, true, false>(
568 Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
569 else
570 eval.template co_gradients<2, false, false, false>(
571 Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
572 }
573 else if constexpr (dim == 2)
574 {
575 if (integration_flag & EvaluationFlags::values)
576 eval.template co_gradients<1, false, true, false>(
577 Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
578 else
579 eval.template co_gradients<1, false, false, false>(
580 Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
581 eval.template co_gradients<0, false, true, false>(
582 Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
583 }
584 else if constexpr (dim == 3)
585 {
586 if (integration_flag & EvaluationFlags::values)
587 eval.template co_gradients<2, false, true, false>(
588 Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
589 else
590 eval.template co_gradients<2, false, false, false>(
591 Kokkos::subview(grad_u, Kokkos::ALL, 2), u);
592 eval.template co_gradients<1, false, true, false>(
593 Kokkos::subview(grad_u, Kokkos::ALL, 1), u);
594 eval.template co_gradients<0, false, true, false>(
595 Kokkos::subview(grad_u, Kokkos::ALL, 0), u);
596 }
597 else
598 Assert(false, ExcMessage("dim must not exceed 3!"));
599 }
600
601 // transform back to the original space
602 if constexpr (dim > 2)
603 eval.template values<2, false, false, true>(u, u);
604 if constexpr (dim > 1)
605 eval.template values<1, false, false, true>(u, u);
606 eval.template values<0, false, false, true>(u, u);
607 }
608 }
609 };
610 } // end of namespace internal
611} // end of namespace Portable
612
613
615
616#endif
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_HOST_DEVICE
Definition config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< index_type > data
Definition mpi.cc:746
EvaluationFlags
The EvaluationFlags enum.
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
constexpr T pow(const T base, const int iexp)
Definition utilities.h:967
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const typename MatrixFree< dim, Number >::Data *data)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const typename MatrixFree< dim, Number >::Data *data)
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const typename MatrixFree< dim, Number >::Data *data)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const typename MatrixFree< dim, Number >::Data *data)
Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type TeamHandle
static void evaluate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags evaluation_flag, const typename MatrixFree< dim, Number >::Data *data)
static void integrate(const unsigned int n_components, const EvaluationFlags::EvaluationFlags integration_flag, const typename MatrixFree< dim, Number >::Data *data)
Kokkos::View< Number *, MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space, Kokkos::MemoryTraits< Kokkos::Unmanaged > > SharedView