deal.II version GIT relicensing-3515-ga7eb338636 2025-06-16 13: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_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
28#include <deal.II/matrix_free/portable_matrix_free.templates.h>
30
31#include <Kokkos_Core.hpp>
32
34
38namespace Portable
39{
63 template <int dim,
64 int fe_degree,
65 int n_q_points_1d = fe_degree + 1,
66 int n_components_ = 1,
67 typename Number = double>
69 {
70 public:
74 using value_type = std::conditional_t<(n_components_ == 1),
75 Number,
77
81 using gradient_type = std::conditional_t<
82 n_components_ == 1,
84 std::conditional_t<n_components_ == dim,
87
92
96 static constexpr unsigned int dimension = dim;
97
101 static constexpr unsigned int n_components = n_components_;
102
106 static constexpr unsigned int n_q_points =
107 Utilities::pow(n_q_points_1d, dim);
108
113 static constexpr unsigned int tensor_dofs_per_component =
114 Utilities::pow(fe_degree + 1, dim);
115
120 static constexpr unsigned int tensor_dofs_per_cell =
122
132 explicit FEEvaluation(const data_type *data,
133 const unsigned int dof_index = 0);
134
139 int
141
148 const data_type *
150
161
170
179 evaluate(const EvaluationFlags::EvaluationFlags evaluate_flag);
180
188 integrate(const EvaluationFlags::EvaluationFlags integration_flag);
189
195 get_value(int q_point) const;
196
202 get_dof_value(int q_point) const;
203
209 submit_value(const value_type &val_in, int q_point);
210
216 submit_dof_value(const value_type &val_in, int q_point);
217
223 get_gradient(int q_point) const;
224
230 submit_gradient(const gradient_type &grad_in, int q_point);
231
232 // clang-format off
243 // clang-format on
244 template <typename Functor>
246 apply_for_each_quad_point(const Functor &func);
247
248 private:
253 };
254
255
256
257 template <int dim,
258 int fe_degree,
259 int n_q_points_1d,
260 int n_components_,
261 typename Number>
264 FEEvaluation(const data_type *data, const unsigned int dof_index)
265 : data(data)
266 , precomputed_data(data->precomputed_data)
267 , shared_data(data->shared_data)
268 , cell_id(data->team_member.league_rank())
269 {
270 AssertIndexRange(dof_index, data->n_dofhandler);
271 }
272
273
274
275 template <int dim,
276 int fe_degree,
277 int n_q_points_1d,
278 int n_components_,
279 typename Number>
286
287
288
289 template <int dim,
290 int fe_degree,
291 int n_q_points_1d,
292 int n_components_,
293 typename Number>
294 DEAL_II_HOST_DEVICE const typename FEEvaluation<dim,
295 fe_degree,
296 n_q_points_1d,
297 n_components_,
298 Number>::data_type *
304
305
306
307 template <int dim,
308 int fe_degree,
309 int n_q_points_1d,
310 int n_components_,
311 typename Number>
315 {
316 // Populate the scratch memory
317 Kokkos::parallel_for(Kokkos::TeamThreadRange(data->team_member,
318 tensor_dofs_per_component),
319 [&](const int &i) {
320 for (unsigned int c = 0; c < n_components_; ++c)
321 shared_data->values(i, c) =
322 src[precomputed_data->local_to_global(
323 i + tensor_dofs_per_component * c, cell_id)];
324 });
325 data->team_member.team_barrier();
326
327 for (unsigned int c = 0; c < n_components_; ++c)
328 {
329 if (precomputed_data->constraint_mask(cell_id * n_components + c) !=
331 unconstrained)
332 internal::resolve_hanging_nodes<dim, fe_degree, false, Number>(
333 data->team_member,
334 precomputed_data->constraint_weights,
335 precomputed_data->constraint_mask(cell_id * n_components + c),
336 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
337 }
338 }
339
340
341
342 template <int dim,
343 int fe_degree,
344 int n_q_points_1d,
345 int n_components_,
346 typename Number>
350 {
351 for (unsigned int c = 0; c < n_components_; ++c)
352 {
353 if (precomputed_data->constraint_mask(cell_id * n_components + c) !=
355 unconstrained)
356 internal::resolve_hanging_nodes<dim, fe_degree, true, Number>(
357 data->team_member,
358 precomputed_data->constraint_weights,
359 precomputed_data->constraint_mask(cell_id * n_components + c),
360 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
361 }
362
363 if (precomputed_data->use_coloring)
364 {
365 Kokkos::parallel_for(
366 Kokkos::TeamThreadRange(data->team_member, tensor_dofs_per_component),
367 [&](const int &i) {
368 for (unsigned int c = 0; c < n_components_; ++c)
369 dst[precomputed_data->local_to_global(
370 i + tensor_dofs_per_component * c, cell_id)] +=
371 shared_data->values(i, c);
372 });
373 }
374 else
375 {
376 Kokkos::parallel_for(
377 Kokkos::TeamThreadRange(data->team_member, tensor_dofs_per_component),
378 [&](const int &i) {
379 for (unsigned int c = 0; c < n_components_; ++c)
380 Kokkos::atomic_add(&dst[precomputed_data->local_to_global(
381 i + (tensor_dofs_per_component)*c, cell_id)],
382 shared_data->values(i, c));
383 });
384 }
385 }
386
387
388
389 template <int dim,
390 int fe_degree,
391 int n_q_points_1d,
392 int n_components,
393 typename Number>
396 const EvaluationFlags::EvaluationFlags evaluation_flag)
397 {
399
400 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
401 precomputed_data->element_type ==
402 ElementType::tensor_symmetric_collocation)
403 {
405 n_components, evaluation_flag, data);
406 }
407 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
408 // shape_info.h for more details
409 else if (fe_degree >= 0 &&
410 internal::use_collocation_evaluation(fe_degree, n_q_points_1d) &&
411 precomputed_data->element_type <= ElementType::tensor_symmetric)
412 {
414 dim,
415 fe_degree,
416 n_q_points_1d,
417 Number>::evaluate(n_components, evaluation_flag, data);
418 }
419 else if (fe_degree >= 0 && precomputed_data->element_type <=
420 ElementType::tensor_symmetric_no_collocation)
421 {
423 evaluate(n_components, evaluation_flag, data);
424 }
425 else
426 {
427 Kokkos::abort("The element type is not yet supported by the portable "
428 "matrix-free module.");
429 }
430 }
431
432
433
434 template <int dim,
435 int fe_degree,
436 int n_q_points_1d,
437 int n_components_,
438 typename Number>
441 const EvaluationFlags::EvaluationFlags integration_flag)
442 {
444
445 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
446 precomputed_data->element_type ==
447 ElementType::tensor_symmetric_collocation)
448 {
450 integrate(n_components, integration_flag, data);
451 }
452 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
453 // shape_info.h for more details
454 else if (fe_degree >= 0 &&
455 internal::use_collocation_evaluation(fe_degree, n_q_points_1d) &&
456 precomputed_data->element_type <= ElementType::tensor_symmetric)
457 {
459 dim,
460 fe_degree,
461 n_q_points_1d,
462 Number>::integrate(n_components, integration_flag, data);
463 }
464 else if (fe_degree >= 0 && precomputed_data->element_type <=
465 ElementType::tensor_symmetric_no_collocation)
466 {
468 integrate(n_components, integration_flag, data);
469 }
470 else
471 {
472 Kokkos::abort("The element type is not yet supported by the portable "
473 "matrix-free module.");
474 }
475 }
476
477
478
479 template <int dim,
480 int fe_degree,
481 int n_q_points_1d,
482 int n_components_,
483 typename Number>
485 fe_degree,
486 n_q_points_1d,
487 n_components_,
488 Number>::value_type
490 int q_point) const
491 {
492 Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
493 if constexpr (n_components_ == 1)
494 {
495 return shared_data->values(q_point, 0);
496 }
497 else
498 {
499 value_type result;
500 for (unsigned int c = 0; c < n_components; ++c)
501 result[c] = shared_data->values(q_point, c);
502 return result;
503 }
504 }
505
506
507
508 template <int dim,
509 int fe_degree,
510 int n_q_points_1d,
511 int n_components_,
512 typename Number>
514 fe_degree,
515 n_q_points_1d,
516 n_components_,
517 Number>::value_type
519 get_dof_value(int dof) const
520 {
521 Assert(dof >= 0 && dof < static_cast<int>(tensor_dofs_per_component),
523 if constexpr (n_components_ == 1)
524 {
525 return shared_data->values(dof, 0);
526 }
527 else
528 {
529 value_type result;
530 for (unsigned int c = 0; c < n_components; ++c)
531 result[c] = shared_data->values(dof, c);
532 return result;
533 }
534 }
535
536
537
538 template <int dim,
539 int fe_degree,
540 int n_q_points_1d,
541 int n_components_,
542 typename Number>
545 submit_value(const value_type &val_in, int q_point)
546 {
547 Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
548 if constexpr (n_components_ == 1)
549 {
550 shared_data->values(q_point, 0) =
551 val_in * precomputed_data->JxW(q_point, cell_id);
552 }
553 else
554 {
555 for (unsigned int c = 0; c < n_components; ++c)
556 shared_data->values(q_point, c) =
557 val_in[c] * precomputed_data->JxW(q_point, cell_id);
558 }
559 }
560
561
562
563 template <int dim,
564 int fe_degree,
565 int n_q_points_1d,
566 int n_components_,
567 typename Number>
570 submit_dof_value(const value_type &val_in, int dof)
571 {
572 Assert(dof >= 0 && dof < tensor_dofs_per_component, ExcInternalError());
573 if constexpr (n_components_ == 1)
574 {
575 shared_data->values(dof, 0) = val_in;
576 }
577 else
578 {
579 for (unsigned int c = 0; c < n_components; ++c)
580 shared_data->values(dof, c) = val_in[c];
581 }
582 }
583
584
585
586 template <int dim,
587 int fe_degree,
588 int n_q_points_1d,
589 int n_components_,
590 typename Number>
592 fe_degree,
593 n_q_points_1d,
594 n_components_,
595 Number>::gradient_type
597 get_gradient(int q_point) const
598 {
599 Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
600 gradient_type grad;
601
602 if constexpr (n_components_ == 1)
603 {
604 for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
605 {
606 Number tmp = 0.;
607 for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
608 tmp +=
609 precomputed_data->inv_jacobian(q_point, cell_id, d_2, d_1) *
610 shared_data->gradients(q_point, d_2, 0);
611 grad[d_1] = tmp;
612 }
613 }
614 else
615 {
616 for (unsigned int c = 0; c < n_components; ++c)
617 for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
618 {
619 Number tmp = 0.;
620 for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
621 tmp +=
622 precomputed_data->inv_jacobian(q_point, cell_id, d_2, d_1) *
623 shared_data->gradients(q_point, d_2, c);
624 grad[c][d_1] = tmp;
625 }
626 }
627
628 return grad;
629 }
630
631
632
633 template <int dim,
634 int fe_degree,
635 int n_q_points_1d,
636 int n_components_,
637 typename Number>
640 submit_gradient(const gradient_type &grad_in, int q_point)
641 {
642 Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
643 if constexpr (n_components_ == 1)
644 {
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 +=
650 precomputed_data->inv_jacobian(q_point, cell_id, d_1, d_2) *
651 grad_in[d_2];
652 shared_data->gradients(q_point, d_1, 0) =
653 tmp * precomputed_data->JxW(q_point, cell_id);
654 }
655 }
656 else
657 {
658 for (unsigned int c = 0; c < n_components; ++c)
659 for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
660 {
661 Number tmp = 0.;
662 for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
663 tmp +=
664 precomputed_data->inv_jacobian(q_point, cell_id, d_1, d_2) *
665 grad_in[c][d_2];
666 shared_data->gradients(q_point, d_1, c) =
667 tmp * precomputed_data->JxW(q_point, cell_id);
668 }
669 }
670 }
671
672
673
674 template <int dim,
675 int fe_degree,
676 int n_q_points_1d,
677 int n_components_,
678 typename Number>
679 template <typename Functor>
682 apply_for_each_quad_point(const Functor &func)
683 {
684 Kokkos::parallel_for(Kokkos::TeamThreadRange(data->team_member, n_q_points),
685 [&](const int &i) { func(this, i); });
686 data->team_member.team_barrier();
687 }
688
689
690
691#ifndef DOXYGEN
692 template <int dim,
693 int fe_degree,
694 int n_q_points_1d,
695 int n_components_,
696 typename Number>
697 constexpr unsigned int
699 n_q_points;
700#endif
701} // namespace Portable
702
704
705#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
const data_type * get_matrix_free_data()
const MatrixFree< dim, Number >::Data * data
static constexpr unsigned int n_components
void evaluate(const EvaluationFlags::EvaluationFlags evaluate_flag)
void read_dof_values(const DeviceVector< Number > &src)
value_type get_value(int q_point) const
void distribute_local_to_global(DeviceVector< Number > &dst) const
static constexpr unsigned int tensor_dofs_per_cell
value_type get_dof_value(int q_point) const
const MatrixFree< dim, Number >::PrecomputedData * precomputed_data
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
gradient_type get_gradient(int q_point) const
static constexpr unsigned int tensor_dofs_per_component
typename MatrixFree< dim, Number >::Data data_type
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:35
#define DEAL_II_HOST_DEVICE
Definition config.h:166
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
std::vector< index_type > data
Definition mpi.cc:750
EvaluationFlags
The EvaluationFlags enum.
constexpr bool use_collocation_evaluation(const unsigned int fe_degree, const unsigned int n_q_points_1d)
Kokkos::View< Number *, MemorySpace::Default::kokkos_space > DeviceVector
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)