Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3083-g7b89508ac7 2025-04-18 12:50: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_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 (void)dof_index;
271 AssertIndexRange(dof_index, data->n_dofhandler);
272 }
273
274
275
276 template <int dim,
277 int fe_degree,
278 int n_q_points_1d,
279 int n_components_,
280 typename Number>
287
288
289
290 template <int dim,
291 int fe_degree,
292 int n_q_points_1d,
293 int n_components_,
294 typename Number>
295 DEAL_II_HOST_DEVICE const typename FEEvaluation<dim,
296 fe_degree,
297 n_q_points_1d,
298 n_components_,
299 Number>::data_type *
305
306
307
308 template <int dim,
309 int fe_degree,
310 int n_q_points_1d,
311 int n_components_,
312 typename Number>
316 {
317 // Populate the scratch memory
318 Kokkos::parallel_for(Kokkos::TeamThreadRange(data->team_member,
319 tensor_dofs_per_component),
320 [&](const int &i) {
321 for (unsigned int c = 0; c < n_components_; ++c)
322 shared_data->values(i, c) =
323 src[precomputed_data->local_to_global(
324 i + tensor_dofs_per_component * c, cell_id)];
325 });
326 data->team_member.team_barrier();
327
328 for (unsigned int c = 0; c < n_components_; ++c)
329 {
330 if (precomputed_data->constraint_mask(cell_id * n_components + c) !=
332 unconstrained)
333 internal::resolve_hanging_nodes<dim, fe_degree, false, Number>(
334 data->team_member,
335 precomputed_data->constraint_weights,
336 precomputed_data->constraint_mask(cell_id * n_components + c),
337 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
338 }
339 }
340
341
342
343 template <int dim,
344 int fe_degree,
345 int n_q_points_1d,
346 int n_components_,
347 typename Number>
351 {
352 for (unsigned int c = 0; c < n_components_; ++c)
353 {
354 if (precomputed_data->constraint_mask(cell_id * n_components + c) !=
356 unconstrained)
357 internal::resolve_hanging_nodes<dim, fe_degree, true, Number>(
358 data->team_member,
359 precomputed_data->constraint_weights,
360 precomputed_data->constraint_mask(cell_id * n_components + c),
361 Kokkos::subview(shared_data->values, Kokkos::ALL, c));
362 }
363
364 if (precomputed_data->use_coloring)
365 {
366 Kokkos::parallel_for(
367 Kokkos::TeamThreadRange(data->team_member, tensor_dofs_per_component),
368 [&](const int &i) {
369 for (unsigned int c = 0; c < n_components_; ++c)
370 dst[precomputed_data->local_to_global(
371 i + tensor_dofs_per_component * c, cell_id)] +=
372 shared_data->values(i, c);
373 });
374 }
375 else
376 {
377 Kokkos::parallel_for(
378 Kokkos::TeamThreadRange(data->team_member, tensor_dofs_per_component),
379 [&](const int &i) {
380 for (unsigned int c = 0; c < n_components_; ++c)
381 Kokkos::atomic_add(&dst[precomputed_data->local_to_global(
382 i + (tensor_dofs_per_component)*c, cell_id)],
383 shared_data->values(i, c));
384 });
385 }
386 }
387
388
389
390 template <int dim,
391 int fe_degree,
392 int n_q_points_1d,
393 int n_components,
394 typename Number>
397 const EvaluationFlags::EvaluationFlags evaluation_flag)
398 {
400
401 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
402 precomputed_data->element_type ==
403 ElementType::tensor_symmetric_collocation)
404 {
406 n_components, evaluation_flag, data);
407 }
408 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
409 // shape_info.h for more details
410 else if (fe_degree >= 0 &&
411 internal::use_collocation_evaluation(fe_degree, n_q_points_1d) &&
412 precomputed_data->element_type <= ElementType::tensor_symmetric)
413 {
415 dim,
416 fe_degree,
417 n_q_points_1d,
418 Number>::evaluate(n_components, evaluation_flag, data);
419 }
420 else if (fe_degree >= 0 && precomputed_data->element_type <=
421 ElementType::tensor_symmetric_no_collocation)
422 {
424 evaluate(n_components, evaluation_flag, data);
425 }
426 else
427 {
428 Kokkos::abort("The element type is not yet supported by the portable "
429 "matrix-free module.");
430 }
431 }
432
433
434
435 template <int dim,
436 int fe_degree,
437 int n_q_points_1d,
438 int n_components_,
439 typename Number>
442 const EvaluationFlags::EvaluationFlags integration_flag)
443 {
445
446 if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d &&
447 precomputed_data->element_type ==
448 ElementType::tensor_symmetric_collocation)
449 {
451 integrate(n_components, integration_flag, data);
452 }
453 // '<=' on type means tensor_symmetric or tensor_symmetric_hermite, see
454 // shape_info.h for more details
455 else if (fe_degree >= 0 &&
456 internal::use_collocation_evaluation(fe_degree, n_q_points_1d) &&
457 precomputed_data->element_type <= ElementType::tensor_symmetric)
458 {
460 dim,
461 fe_degree,
462 n_q_points_1d,
463 Number>::integrate(n_components, integration_flag, data);
464 }
465 else if (fe_degree >= 0 && precomputed_data->element_type <=
466 ElementType::tensor_symmetric_no_collocation)
467 {
469 integrate(n_components, integration_flag, data);
470 }
471 else
472 {
473 Kokkos::abort("The element type is not yet supported by the portable "
474 "matrix-free module.");
475 }
476 }
477
478
479
480 template <int dim,
481 int fe_degree,
482 int n_q_points_1d,
483 int n_components_,
484 typename Number>
486 fe_degree,
487 n_q_points_1d,
488 n_components_,
489 Number>::value_type
491 int q_point) const
492 {
493 Assert(q_point >= 0 && q_point < n_q_points, ExcInternalError());
494 if constexpr (n_components_ == 1)
495 {
496 return shared_data->values(q_point, 0);
497 }
498 else
499 {
500 value_type result;
501 for (unsigned int c = 0; c < n_components; ++c)
502 result[c] = shared_data->values(q_point, c);
503 return result;
504 }
505 }
506
507
508
509 template <int dim,
510 int fe_degree,
511 int n_q_points_1d,
512 int n_components_,
513 typename Number>
515 fe_degree,
516 n_q_points_1d,
517 n_components_,
518 Number>::value_type
520 get_dof_value(int dof) const
521 {
522 Assert(dof >= 0 && dof < tensor_dofs_per_component, ExcInternalError());
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:746
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)