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, typename Val>
      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, Val val) {
      51              :           // The number of indices is equal to the view rank
      52              : 
      53              :           /** If scatterToPoint ultimately writes into `view`, convert once to its element type.
      54              :               Now Val can c.f be double or floar
      55              :           */
      56              :           using out_t =
      57              :             std::remove_cv_t<std::remove_reference_t<typename View::non_const_value_type>>;
      58         4992 :           const out_t v_out = static_cast<out_t>(val);
      59         4992 :             (scatterToPoint<ScatterPoint>(std::make_index_sequence<View::rank>{}, view, wlo, whi,
      60              :                                           args, v_out),
      61              :              ...);
      62         4992 :         }
      63              : 
      64              :         template <unsigned long GatherPoint, unsigned long... Index, typename View, typename T,
      65              :                   typename IndexType>
      66         8064 :         KOKKOS_INLINE_FUNCTION constexpr typename View::value_type gatherFromPoint(
      67              :             const std::index_sequence<Index...>&, const View& view,
      68              :             const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
      69              :             const Vector<IndexType, View::rank>& args) {
      70         8064 :             return (interpolationWeight<GatherPoint, Index>(wlo, whi) * ...)
      71        16128 :                    * view(interpolationIndex<GatherPoint, Index>(args)...);
      72              :         }
      73              : 
      74              :         template <unsigned long... GatherPoint, typename View, typename T, typename IndexType>
      75          384 :         KOKKOS_INLINE_FUNCTION constexpr typename View::value_type gatherFromField(
      76              :             const std::index_sequence<GatherPoint...>&, const View& view,
      77              :             const Vector<T, View::rank>& wlo, const Vector<T, View::rank>& whi,
      78              :             const Vector<IndexType, View::rank>& args) {
      79              :             // The number of indices is equal to the view rank
      80          384 :             return (gatherFromPoint<GatherPoint>(std::make_index_sequence<View::rank>{}, view, wlo,
      81              :                                                  whi, args)
      82          384 :                     + ...);
      83              :         }
      84              :     }  // namespace detail
      85              : }  // namespace ippl