BinningTools.h

Go to the documentation of this file.

#ifndef BINNINGTOOLS_H
#define BINNINGTOOLS_H
 
#include "ParallelReduceTools.h"  // needed for HistoReductionMode and maxArrSize
#include "Utilities/OpalException.h"
 
namespace ParticleBinning {
 
    template <typename bin_index_type>
    HistoReductionMode determineHistoReductionMode(
            HistoReductionMode modePreference, bin_index_type binCount) {
        // Overwrite standard mode if compiled with default host execution space!
        if (std::is_same<Kokkos::DefaultExecutionSpace, Kokkos::DefaultHostExecutionSpace>::value)
            return HistoReductionMode::HostOnly;
 
#ifdef OPALX_DEVICE_COMPILATION
        // In device builds HostOnly reductions are not supported; calling them is a logic error.
        if (modePreference == HistoReductionMode::HostOnly) {
            throw OpalException(
                    "ParticleBinning::determineHistoReductionMode",
                    "HistoReductionMode::HostOnly is not supported when OPALX_DEVICE_COMPILATION "
                    "is "
                    "enabled.");
        }
#endif
        // Otherwise choose automatically if Standard and respect preference if not on host and not
        // standard!
        if (modePreference == HistoReductionMode::Standard) {
            return (binCount <= maxArrSize<bin_index_type>) ? HistoReductionMode::ParallelReduce
                                                            : HistoReductionMode::TeamBased;
        } else {
            return modePreference;
        }
    }
 
    template <typename bunch_type>
    struct CoordinateSelector {
        using value_type = typename bunch_type::Layout_t::value_type;
 
        using size_type = typename bunch_type::size_type;
 
        using position_view_type = typename bunch_type::particle_position_type::view_type;
 
        position_view_type data_arr;  
        const int axis;               
 
        CoordinateSelector(int axis_) : axis(axis_) {}
 
        void updateDataArr(bunch_type& bunch) { data_arr = bunch.P.getView(); }
 
        KOKKOS_INLINE_FUNCTION
        value_type operator()(const size_type& i) const {
            // std::cout << "CoordinateSelector: " << i << std::endl; // TODO: debug, remove later
            const value_type value = fabs(data_arr(i)[axis]);
            return value / sqrt(1 + value * value);  // Normalize to v/c, so v in [0, 1]
        }
    };
 
    template <typename bunch_type>
    struct GammaSelector {
        using value_type         = typename bunch_type::Layout_t::value_type;
        using size_type          = typename bunch_type::size_type;
        using position_view_type = typename bunch_type::particle_position_type::view_type;
 
        position_view_type data_arr;
        const int axis;
 
        GammaSelector(int axis_) : axis(axis_) {}
 
        void updateDataArr(bunch_type& bunch) { data_arr = bunch.P.getView(); }
 
        KOKKOS_INLINE_FUNCTION
        value_type operator()(const size_type& i) const {
            const value_type value = data_arr(i)[axis];
            return sqrt(1 + value * value);
        }
    };
 
    template <typename ViewType>
    void computeFixSum(const ViewType& input_view, const ViewType& post_sum_view) {
        using execution_space = typename ViewType::execution_space;
        using size_type       = typename ViewType::size_type;
        using value_type      = typename ViewType::value_type;
 
        // Ensure the output view has the correct size
        if (post_sum_view.extent(0) != input_view.extent(0) + 1) {
            Inform m("computePostSum");
            m << level4 << "Output view must have size input_view.extent(0) + 1" << endl;
            ippl::Comm->abort();
        }
 
        // Initialize the first element to 0
        Kokkos::parallel_for(
                "InitPostSum", Kokkos::RangePolicy<execution_space>(0, 1),
                KOKKOS_LAMBDA(const size_type) { post_sum_view(0) = 0; });
 
        // Compute the fix sum
        Kokkos::parallel_scan(
                "ComputePostSum", Kokkos::RangePolicy<execution_space>(0, input_view.extent(0)),
                KOKKOS_LAMBDA(const size_type& i, value_type& partial_sum, bool final) {
                    partial_sum += input_view(i);
                    if (final) {
                        post_sum_view(i + 1) = partial_sum;
                    }
                });
    }
 
    template <typename ValueType, typename SizeType, typename HashType>
    bool viewIsSorted(const Kokkos::View<ValueType*> view, HashType indices, SizeType npart) {
        bool sorted = true;
        Kokkos::parallel_reduce(
                "CheckSorted", npart - 1,
                KOKKOS_LAMBDA(const SizeType& i, bool& update) {
                    if (view(indices(i)) > view(indices(i + 1))) update = false;
                },
                Kokkos::LAnd<bool>(sorted));
        return sorted;
    }
 
}  // namespace ParticleBinning
 
#endif  // BINNINGTOOLS_H