Line data Source code
1 : //
2 : // Class CIC
3 : // First order/cloud-in-cell grid interpolation. Currently implemented as
4 : // global functions, but in order to support higher or lower order interpolation,
5 : // these should be moved into structs.
6 : //
7 :
8 : namespace ippl {
9 : namespace detail {
10 : template <unsigned long Point, unsigned long Index, typename Weights>
11 575232 : KOKKOS_INLINE_FUNCTION constexpr typename Weights::value_type interpolationWeight(
12 : const Weights& wlo, const Weights& whi) {
13 : if constexpr (Point & (1 << Index)) {
14 287616 : return wlo[Index];
15 : } else {
16 287616 : return whi[Index];
17 : }
18 : // device code cannot throw exceptions, but we need a
19 : // dummy return to silence the warning
20 : return 0;
21 : }
22 :
23 : template <unsigned long Point, unsigned long Index, typename Indices>
24 575232 : KOKKOS_INLINE_FUNCTION constexpr typename Indices::value_type interpolationIndex(
25 : const Indices& args) {
26 : if constexpr (Point & (1 << Index)) {
27 287616 : return args[Index] - 1;
28 : } else {
29 287616 : return args[Index];
30 : }
31 : // device code cannot throw exceptions, but we need a
32 : // dummy return to silence the warning
33 : return 0;
34 : }
35 :
36 : template <unsigned long ScatterPoint, unsigned long... Index, typename View, typename T,
37 : typename IndexType>
38 104832 : KOKKOS_INLINE_FUNCTION constexpr void scatterToPoint(
39 : const std::index_sequence<Index...>&, const View& view,
40 : const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
41 : const Vector<IndexType, View::rank>& args, const T& val) {
42 209664 : Kokkos::atomic_add(&view(interpolationIndex<ScatterPoint, Index>(args)...),
43 104832 : val * (interpolationWeight<ScatterPoint, Index>(wlo, whi) * ...));
44 104832 : }
45 :
46 : template <unsigned long... ScatterPoint, typename View, typename T, typename IndexType>
47 4992 : KOKKOS_INLINE_FUNCTION constexpr void scatterToField(
48 : const std::index_sequence<ScatterPoint...>&, const View& view,
49 : const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
50 : const Vector<IndexType, View::rank>& args, T val) {
51 : // The number of indices is equal to the view rank
52 4992 : (scatterToPoint<ScatterPoint>(std::make_index_sequence<View::rank>{}, view, wlo, whi,
53 : args, val),
54 : ...);
55 4992 : }
56 :
57 : template <unsigned long GatherPoint, unsigned long... Index, typename View, typename T,
58 : typename IndexType>
59 8064 : KOKKOS_INLINE_FUNCTION constexpr typename View::value_type gatherFromPoint(
60 : const std::index_sequence<Index...>&, const View& view,
61 : const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
62 : const Vector<IndexType, View::rank>& args) {
63 8064 : return (interpolationWeight<GatherPoint, Index>(wlo, whi) * ...)
64 16128 : * view(interpolationIndex<GatherPoint, Index>(args)...);
65 : }
66 :
67 : template <unsigned long... GatherPoint, typename View, typename T, typename IndexType>
68 384 : KOKKOS_INLINE_FUNCTION constexpr typename View::value_type gatherFromField(
69 : const std::index_sequence<GatherPoint...>&, const View& view,
70 : const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
71 : const Vector<IndexType, View::rank>& args) {
72 : // The number of indices is equal to the view rank
73 384 : return (gatherFromPoint<GatherPoint>(std::make_index_sequence<View::rank>{}, view, wlo,
74 : whi, args)
75 384 : + ...);
76 : }
77 : } // namespace detail
78 : } // namespace ippl
|