deal.II version GIT relicensing-2594-g94a7da5c31 2025-02-10 10:20: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_fe_evaluation.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2023 - 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_fe_evaluation_h
16#define dealii_portable_fe_evaluation_h
17
18#include <deal.II/base/config.h>
19
21#include <deal.II/base/tensor.h>
23
27#include <deal.II/matrix_free/portable_matrix_free.templates.h>
29
30#include <Kokkos_Core.hpp>
31
33
37namespace Portable
38{
62 template <int dim,
63 int fe_degree,
64 int n_q_points_1d = fe_degree + 1,
65 int n_components_ = 1,
66 typename Number = double>
68 {
69 public:
73 using value_type = std::conditional_t<(n_components_ == 1),
74 Number,
76
80 using gradient_type = std::conditional_t<
81 n_components_ == 1,
83 std::conditional_t<n_components_ == dim,
86
91
95 static constexpr unsigned int dimension = dim;
96
100 static constexpr unsigned int n_components = n_components_;
101
105 static constexpr unsigned int n_q_points =
106 Utilities::pow(n_q_points_1d, dim);
107
111 static constexpr unsigned int tensor_dofs_per_cell =
112 Utilities::pow(fe_degree + 1, dim) * n_components;
113
119
129 read_dof_values(const Number *src);
130
138 distribute_local_to_global(Number *dst) const;
139
148 evaluate(const EvaluationFlags::EvaluationFlags evaluate_flag);
149
157 DEAL_II_DEPRECATED_WITH_COMMENT("Use the version taking EvaluationFlags.")
159 void
160 evaluate(const bool evaluate_val, const bool evaluate_grad);
161
169 integrate(const EvaluationFlags::EvaluationFlags integration_flag);
170
178 DEAL_II_DEPRECATED_WITH_COMMENT("Use the version taking EvaluationFlags.")
180 void
181 integrate(const bool integrate_val, const bool integrate_grad);
182
188 get_value(int q_point) const;
189
195 get_dof_value(int q_point) const;
196
202 submit_value(const value_type &val_in, int q_point);
203
209 submit_dof_value(const value_type &val_in, int q_point);
210
216 get_gradient(int q_point) const;
217
223 submit_gradient(const gradient_type &grad_in, int q_point);
224
225 // clang-format off
236 // clang-format on
237 template <typename Functor>
239 apply_for_each_quad_point(const Functor &func);
240
241 private:
243 SharedData<dim, Number> *shared_data;
245 };
246
247
248
249 template <int dim,
250 int fe_degree,
251 int n_q_points_1d,
252 int n_components_,
253 typename Number>
255 FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
256 FEEvaluation(const data_type *data, SharedData<dim, Number> *shdata)
257 : data(data)
258 , shared_data(shdata)
259 , cell_id(shared_data->team_member.league_rank())
260 {}
261
262
263
264 template <int dim,
265 int fe_degree,
266 int n_q_points_1d,
267 int n_components_,
268 typename Number>
271 read_dof_values(const Number *src)
272 {
273 // Populate the scratch memory
274 Kokkos::parallel_for(
275 Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
276 [&](const int &i) {
277 for (unsigned int c = 0; c < n_components_; ++c)
278 shared_data->values(i, c) = src[data->local_to_global(
279 cell_id, i + (tensor_dofs_per_cell / n_components) * c)];
280 });
281 shared_data->team_member.team_barrier();
282
283 for (unsigned int c = 0; c < n_components_; ++c)
284 {
285 internal::resolve_hanging_nodes<dim, fe_degree, false, Number>(
286 shared_data->team_member,
287 data->constraint_weights,
288 data->constraint_mask(cell_id * n_components + c),
289 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
290 }
291 }
292
293
294
295 template <int dim,
296 int fe_degree,
297 int n_q_points_1d,
298 int n_components_,
299 typename Number>
302 distribute_local_to_global(Number *dst) const
303 {
304 for (unsigned int c = 0; c < n_components_; ++c)
305 {
306 internal::resolve_hanging_nodes<dim, fe_degree, true, Number>(
307 shared_data->team_member,
308 data->constraint_weights,
309 data->constraint_mask(cell_id * n_components + c),
310 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
311 }
312
313 if (data->use_coloring)
314 {
315 Kokkos::parallel_for(
316 Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
317 [&](const int &i) {
318 for (unsigned int c = 0; c < n_components_; ++c)
319 dst[data->local_to_global(
320 cell_id, i + (tensor_dofs_per_cell / n_components) * c)] +=
321 shared_data->values(i, c);
322 });
323 }
324 else
325 {
326 Kokkos::parallel_for(
327 Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
328 [&](const int &i) {
329 for (unsigned int c = 0; c < n_components_; ++c)
330 Kokkos::atomic_add(
331 &dst[data->local_to_global(
332 cell_id, i + (tensor_dofs_per_cell / n_components) * c)],
333 shared_data->values(i, c));
334 });
335 }
336 }
337
338
339
340 template <int dim,
341 int fe_degree,
342 int n_q_points_1d,
343 int n_components_,
344 typename Number>
347 const EvaluationFlags::EvaluationFlags evaluate_flag)
348 {
349 // First evaluate the gradients because it requires values that will be
350 // changed if evaluate_val is true
353 dim,
354 fe_degree,
355 n_q_points_1d,
356 Number>
357 evaluator_tensor_product(shared_data->team_member,
358 data->shape_values,
359 data->shape_gradients,
360 data->co_shape_gradients);
361
362 for (unsigned int c = 0; c < n_components_; ++c)
363 {
364 if ((evaluate_flag & EvaluationFlags::values) &&
365 (evaluate_flag & EvaluationFlags::gradients))
366 {
367 evaluator_tensor_product.evaluate_values_and_gradients(
368 Kokkos::subview(shared_data->values, Kokkos::ALL, c),
369 Kokkos::subview(
370 shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
371 shared_data->team_member.team_barrier();
372 }
373 else if (evaluate_flag & EvaluationFlags::gradients)
374 {
375 evaluator_tensor_product.evaluate_gradients(
376 Kokkos::subview(shared_data->values, Kokkos::ALL, c),
377 Kokkos::subview(
378 shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
379 shared_data->team_member.team_barrier();
380 }
381 else if (evaluate_flag & EvaluationFlags::values)
382 {
383 evaluator_tensor_product.evaluate_values(
384 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
385 shared_data->team_member.team_barrier();
386 }
387 }
388 }
389
390
391
392 template <int dim,
393 int fe_degree,
394 int n_q_points_1d,
395 int n_components_,
396 typename Number>
399 const bool evaluate_val,
400 const bool evaluate_grad)
401 {
402 evaluate(
405 }
406
407
408
409 template <int dim,
410 int fe_degree,
411 int n_q_points_1d,
412 int n_components_,
413 typename Number>
416 const EvaluationFlags::EvaluationFlags integration_flag)
417 {
420 dim,
421 fe_degree,
422 n_q_points_1d,
423 Number>
424 evaluator_tensor_product(shared_data->team_member,
425 data->shape_values,
426 data->shape_gradients,
427 data->co_shape_gradients);
428
429
430 for (unsigned int c = 0; c < n_components_; ++c)
431 {
432 if ((integration_flag & EvaluationFlags::values) &&
433 (integration_flag & EvaluationFlags::gradients))
434 {
435 evaluator_tensor_product.integrate_values_and_gradients(
436 Kokkos::subview(shared_data->values, Kokkos::ALL, c),
437 Kokkos::subview(
438 shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
439 }
440 else if (integration_flag & EvaluationFlags::values)
441 {
442 evaluator_tensor_product.integrate_values(
443 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
444 shared_data->team_member.team_barrier();
445 }
446 else if (integration_flag & EvaluationFlags::gradients)
447 {
448 evaluator_tensor_product.template integrate_gradients<false>(
449 Kokkos::subview(shared_data->values, Kokkos::ALL, c),
450 Kokkos::subview(
451 shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
452 shared_data->team_member.team_barrier();
453 }
454 }
455 }
456
457
458
459 template <int dim,
460 int fe_degree,
461 int n_q_points_1d,
462 int n_components_,
463 typename Number>
466 const bool integrate_val,
467 const bool integrate_grad)
468 {
469 integrate(
472 }
473
474
475
476 template <int dim,
477 int fe_degree,
478 int n_q_points_1d,
479 int n_components_,
480 typename Number>
482 fe_degree,
483 n_q_points_1d,
484 n_components_,
485 Number>::value_type
487 int q_point) const
488 {
489 if constexpr (n_components_ == 1)
490 {
491 return shared_data->values(q_point, 0);
492 }
493 else
494 {
495 value_type result;
496 for (unsigned int c = 0; c < n_components; ++c)
497 result[c] = shared_data->values(q_point, c);
498 return result;
499 }
500 }
501
502
503
504 template <int dim,
505 int fe_degree,
506 int n_q_points_1d,
507 int n_components_,
508 typename Number>
510 fe_degree,
511 n_q_points_1d,
512 n_components_,
513 Number>::value_type
515 get_dof_value(int q_point) const
516 {
517 if constexpr (n_components_ == 1)
518 {
519 return shared_data->values(q_point, 0);
520 }
521 else
522 {
523 value_type result;
524 for (unsigned int c = 0; c < n_components; ++c)
525 result[c] = shared_data->values(q_point, c);
526 return result;
527 }
528 }
529
530
531
532 template <int dim,
533 int fe_degree,
534 int n_q_points_1d,
535 int n_components_,
536 typename Number>
539 submit_value(const value_type &val_in, int q_point)
540 {
541 if constexpr (n_components_ == 1)
542 {
543 shared_data->values(q_point, 0) = val_in * data->JxW(cell_id, q_point);
544 }
545 else
546 {
547 for (unsigned int c = 0; c < n_components; ++c)
548 shared_data->values(q_point, c) =
549 val_in[c] * data->JxW(cell_id, q_point);
550 }
551 }
552
553
554
555 template <int dim,
556 int fe_degree,
557 int n_q_points_1d,
558 int n_components_,
559 typename Number>
562 submit_dof_value(const value_type &val_in, int q_point)
563 {
564 if constexpr (n_components_ == 1)
565 {
566 shared_data->values(q_point, 0) = val_in;
567 }
568 else
569 {
570 for (unsigned int c = 0; c < n_components; ++c)
571 shared_data->values(q_point, c) = val_in[c];
572 }
573 }
574
575
576
577 template <int dim,
578 int fe_degree,
579 int n_q_points_1d,
580 int n_components_,
581 typename Number>
583 fe_degree,
584 n_q_points_1d,
585 n_components_,
586 Number>::gradient_type
588 get_gradient(int q_point) const
589 {
590 gradient_type grad;
591
592 if constexpr (n_components_ == 1)
593 {
594 for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
595 {
596 Number tmp = 0.;
597 for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
598 tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) *
599 shared_data->gradients(q_point, d_2, 0);
600 grad[d_1] = tmp;
601 }
602 }
603 else
604 {
605 for (unsigned int c = 0; c < n_components; ++c)
606 for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
607 {
608 Number tmp = 0.;
609 for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
610 tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) *
611 shared_data->gradients(q_point, d_2, c);
612 grad[c][d_1] = tmp;
613 }
614 }
615
616 return grad;
617 }
618
619
620
621 template <int dim,
622 int fe_degree,
623 int n_q_points_1d,
624 int n_components_,
625 typename Number>
628 submit_gradient(const gradient_type &grad_in, int q_point)
629 {
630 if constexpr (n_components_ == 1)
631 {
632 for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
633 {
634 Number tmp = 0.;
635 for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
636 tmp +=
637 data->inv_jacobian(cell_id, q_point, d_1, d_2) * grad_in[d_2];
638 shared_data->gradients(q_point, d_1, 0) =
639 tmp * data->JxW(cell_id, q_point);
640 }
641 }
642 else
643 {
644 for (unsigned int c = 0; c < n_components; ++c)
645 for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
646 {
647 Number tmp = 0.;
648 for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
649 tmp += data->inv_jacobian(cell_id, q_point, d_1, d_2) *
650 grad_in[c][d_2];
651 shared_data->gradients(q_point, d_1, c) =
652 tmp * data->JxW(cell_id, q_point);
653 }
654 }
655 }
656
657
658
659 template <int dim,
660 int fe_degree,
661 int n_q_points_1d,
662 int n_components_,
663 typename Number>
664 template <typename Functor>
667 apply_for_each_quad_point(const Functor &func)
668 {
669 Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
670 n_q_points),
671 [&](const int &i) { func(this, i); });
672 shared_data->team_member.team_barrier();
673 }
674
675
676
677#ifndef DOXYGEN
678 template <int dim,
679 int fe_degree,
680 int n_q_points_1d,
681 int n_components_,
682 typename Number>
683 constexpr unsigned int
686#endif
687} // namespace Portable
688
690
691#endif
static constexpr unsigned int n_q_points
SharedData< dim, Number > * shared_data
void integrate(const EvaluationFlags::EvaluationFlags integration_flag)
static constexpr unsigned int dimension
static constexpr unsigned int n_components
value_type get_value(int q_point) const
static constexpr unsigned int tensor_dofs_per_cell
value_type get_dof_value(int q_point) const
std::conditional_t< n_components_==1, Tensor< 1, dim, Number >, std::conditional_t< n_components_==dim, Tensor< 2, dim, Number >, Tensor< 1, n_components_, Tensor< 1, dim, Number > > > > gradient_type
std::conditional_t<(n_components_==1), Number, Tensor< 1, n_components_, Number > > value_type
void evaluate(const EvaluationFlags::EvaluationFlags evaluate_flag)
void distribute_local_to_global(Number *dst) const
gradient_type get_gradient(int q_point) const
typename MatrixFree< dim, Number >::Data data_type
void read_dof_values(const Number *src)
void submit_dof_value(const value_type &val_in, int q_point)
void submit_value(const value_type &val_in, int q_point)
void submit_gradient(const gradient_type &grad_in, int q_point)
void apply_for_each_quad_point(const Functor &func)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_DEPRECATED_WITH_COMMENT(comment)
Definition config.h:229
#define DEAL_II_HOST_DEVICE
Definition config.h:114
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
The namespace for the EvaluationFlags enum.
EvaluationFlags
The EvaluationFlags enum.
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966