Add fusion implementation similar to: doi.org/10.1063/5.0051178

Author: FilipO28555

include/picongpu/param/fusion.param

@@ -0,0 +1,103 @@
+/* Copyright 2025 Filip Optolowicz 
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "picongpu/param/particleFilters.param"
+#include "picongpu/particles/fusion/fusion.def"
+#include "picongpu/plugins/misc/SpeciesFilter.hpp"
+
+#include <pmacc/meta/conversion/MakeSeq.hpp>
+
+
+#ifndef PARAM_DYNAMIC_LOG
+#    define PARAM_DYNAMIC_LOG 0
+#endif
+
+namespace picongpu
+{
+    namespace particles
+    {
+        namespace fusion
+        {
+            namespace precision
+            {
+                using float_COLL = float_64;
+            } // namespace precision
+
+            using namespace picongpu::plugins;
+
+            // Example of defining a fusion reaction
+            // Deuterium-Tritium reaction: D + T -> n + He4
+            struct DT
+            {
+
+                using reactants = pmacc::mp_list<
+                    Pair<PIC_Tritons, PIC_Deuterons>
+                    >;
+
+                using products = pmacc::mp_list<
+                    Pair<PIC_Neutrons, PIC_He4>
+                    >;
+
+                struct Params
+                {
+                    // B_g - Gamov constant: B_g = pi*alpha*Z1*Z2*sqrt(2*m_r*c^2) in keV^1/2
+                    static constexpr float_X BG = 34.3827_X;
+
+                    // Coefficients for the Bosch-Hale parameterization of the DT fusion cross-section
+                    static constexpr float_X A1 = 6.927e4_X;
+                    static constexpr float_X A2 = 7.454e8_X;
+                    static constexpr float_X A3 = 2.050e6_X;
+                    static constexpr float_X A4 = 5.2002e4_X;
+                    static constexpr float_X A5 = 0.0_X;
+
+                    static constexpr float_X B1 = 6.38e1_X;
+                    static constexpr float_X B2 = -9.95e-1_X;
+                    static constexpr float_X B3 = 6.981e-5_X;
+                    static constexpr float_X B4 = 1.728e-4_X;
+                };
+
+                using CrossSectionInterpolator = relativistic::FusionFunctor<Params>;
+            };
+
+
+            /** FusionPipeline defines in which order species interact with each other
+             *
+             * the functors are called in order (from first to last functor)
+             */
+            using FusionPipeline = pmacc::mp_list<
+                Collider<DT::CrossSectionInterpolator, DT::reactants, DT::products, OneFilter<filter::All>>>;
+
+            // using FusionPipeline = pmacc::mp_list<ColliderFromStruct<DT>>;
+
+            /** Chunk size used for cell list allocations.
+             *
+             * To reduce the fragmentation of the heap memory on accelerators the collision algorithm is allocating a
+             * multiple of this value to store a cell list of particle IDs. The value must be non zero.
+             */
+            constexpr uint32_t cellListChunkSize = TYPICAL_PARTICLES_PER_CELL;
+            // if weighting of the products were to be less or equal than productMinWeighting we set  Fmult = std::max(1,productWeighting/productMinWeighting);
+            constexpr float_X productMinWeighting = 16.1;
+            constexpr uint32_t maxFmult = 1e6;
+            constexpr bool debugFusion = false;
+            constexpr bool alwaysFuseQ = false;
+        } // namespace fusion
+    } // namespace particles
+} // namespace picongpu

include/picongpu/particles/fusion/Collider.def

@@ -0,0 +1,82 @@
+/* Copyright 2019-2024 Rene Widera, Pawel Ordyna
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "picongpu/particles/filter/All.hpp"
+
+#include <pmacc/meta/Pair.hpp>
+
+namespace picongpu::particles::fusion
+{
+    /* A pair of colliding species.
+     *
+     * @tparam T_Species0
+     * @tparam T_Species1
+     */
+    template<typename T_Species0, typename T_Species1>
+    using Pair = pmacc::meta::Pair<T_Species0, T_Species1>;
+
+    /* A set of particle filters for a pair of colliding species.
+     *
+     * @tparam T_Filter0 Filter applied on the 1st species in a pair.
+     * @tparam T_Filter1 Filter applied on the 2nd species in a pair.
+     */
+    template<typename T_Filter0, typename T_Filter1>
+    using FilterPair = pmacc::meta::Pair<T_Filter0, T_Filter1>;
+
+    /* Sets the same particle filter for both species in a pair.
+     *
+     * @tparam T_Filter A common particle filter for both colliding species.
+     */
+    template<typename T_Filter>
+    using OneFilter = pmacc::meta::Pair<T_Filter, T_Filter>;
+
+    /* Defines a set of binary collisions with common parameters.
+     *
+     * @tparam T_CrossSectionInterpolator A binary particle functor defining a single
+     *     macro particle collision in the binary-collision algorithm.
+     * @tparam T_SpeciesPairs A sequence of pairs of colliding species.
+     * @tparam T_Params A struct defining `coulombLog` for the collisions.
+     * @tparam T_FilterPair A pair of particle filters. Each for every species
+     *     in a pair of colliding species. This pair of filters will be
+     *     aplied to all pairs in T_SpeciesPairs.
+     */
+    template<
+        typename T_CrossSectionInterpolator,
+        typename T_SpeciesPairsReactants,
+        typename T_SpeciesPairsProducts,
+        typename T_FilterPair = OneFilter<filter::All>>
+    struct Collider
+    {
+        using Functor = T_CrossSectionInterpolator;
+        using SpeciesPairsReactants = T_SpeciesPairsReactants;
+        using SpeciesPairsProducts = T_SpeciesPairsProducts;
+        using FilterPair = T_FilterPair;
+    };
+
+    template<typename T_FusionStruct>
+    struct ColliderFromStruct
+    {
+        using Functor = typename T_FusionStruct::CrossSectionInterpolator;
+        using SpeciesPairsReactants = typename T_FusionStruct::reactants;
+        using SpeciesPairsProducts = typename T_FusionStruct::products;
+        using FilterPair = typename T_FusionStruct::FilterPair;
+    };
+} // namespace picongpu::particles::fusion

include/picongpu/particles/fusion/Collider.hpp

@@ -0,0 +1,79 @@
+/* Copyright 2019-2024 Rene Widera, Pawel Ordyna, Filip Optolowicz
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include "picongpu/defines.hpp"
+#include "picongpu/particles/fusion/Collider.def"
+#include "picongpu/particles/fusion/WithPeer.hpp"
+
+#include <pmacc/meta/ForEach.hpp>
+#include <pmacc/meta/accessors/First.hpp>
+#include <pmacc/meta/accessors/Second.hpp>
+#include <pmacc/meta/conversion/ApplyGuard.hpp>
+#include <pmacc/meta/conversion/ToSeq.hpp>
+
+namespace picongpu::particles::fusion
+{
+    namespace detail
+    {
+        // "For each" implementation for calling a collider for each species pair in a list of reactants
+        template<
+            typename T_SpeciesPairListReactants,
+            typename T_SpeciesPairListProducts,
+            typename T_Collider,
+            uint32_t colliderId>
+        struct CallColliderForAPair
+        {
+            template<size_t... I>
+            HINLINE void operator()(
+                std::index_sequence<I...>,
+                std::shared_ptr<DeviceHeap> const& deviceHeap,
+                uint32_t currentStep)
+            {
+                (fusion::WithPeer<
+                     typename T_Collider::Functor,
+                     typename pmacc::mp_at_c<T_SpeciesPairListReactants, I>::first,
+                     typename pmacc::mp_at_c<T_SpeciesPairListReactants, I>::second,
+                     typename pmacc::mp_at_c<T_SpeciesPairListProducts, I>::first,
+                     typename pmacc::mp_at_c<T_SpeciesPairListProducts, I>::second,
+                     typename T_Collider::FilterPair,
+                     colliderId,
+                     I>{}(deviceHeap, currentStep),
+                 ...);
+            }
+        };
+    } // namespace detail
+
+    template<typename T_Collider, uint32_t colliderId>
+    struct CallCollider
+    {
+        void operator()(std::shared_ptr<DeviceHeap> const& deviceHeap, uint32_t currentStep)
+        {
+            using SpeciesPairListReactants = pmacc::ToSeq<typename T_Collider::SpeciesPairsReactants>;
+            using SpeciesPairListProducts = pmacc::ToSeq<typename T_Collider::SpeciesPairsProducts>;
+            constexpr size_t numPairs = pmacc::mp_size<SpeciesPairListReactants>::value;
+            std::make_index_sequence<numPairs> index{};
+            detail::CallColliderForAPair<SpeciesPairListReactants, SpeciesPairListProducts, T_Collider, colliderId>{}(
+                index,
+                deviceHeap,
+                currentStep);
+        }
+    };
+} // namespace picongpu::particles::fusion

include/picongpu/particles/fusion/IBinary.def

@@ -0,0 +1,48 @@
+/* Copyright 2019-2024 Rene Widera, Pawel Ordyna
+ *
+ * This file is part of PIConGPU.
+ *
+ * PIConGPU is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PIConGPU is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with PIConGPU.
+ * If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <pmacc/functor/Interface.hpp>
+
+namespace picongpu
+{
+    namespace particles
+    {
+        namespace fusion
+        {
+            namespace acc
+            {
+                template<typename T_BinaryFunctor>
+                using IBinary = pmacc::functor::acc::Interface<T_BinaryFunctor, 2u, void>;
+            } // namespace acc
+
+            /**  interface for a binary particle functor
+             *
+             *
+             * @tparam T_BinaryFunctor binary particle acc functor, must contain
+             *         `void operator()(P1 & particle1, P2 & particle2, ...)`
+             *          and support at least two particles
+             */
+            template<typename T_BinaryFunctor>
+            using IBinary = pmacc::functor::Interface<T_BinaryFunctor, 2u, void>;
+
+        } // namespace fusion
+    } // namespace particles
+} // namespace picongpu