instant field ionization transitions

include/picongpu/param/atomicPhysics.param

@@ -32,27 +32,38 @@
 
 namespace picongpu::atomicPhysics
 {
-    // configuration of electron ElectronHistogram
-    //{
+    //! configuration of electron ElectronHistogram
+    //!@{
     struct MaxEnergyParam
     {
-        // eV
+        // unit: eV
         static constexpr float_X maxEnergy = float_X(1e7);
     };
     using ElectronHistogram = particles::atomicPhysics::electronDistribution::LogSpaceHistogram<
-        100u, // T_numberBins
+        // T_numberBins
+        100u,
         MaxEnergyParam>;
-    //}
+    //!@}
 
     //! atomicPhysics rate solver settings
     struct RateSolverParam
     {
-        // atomicPhysics timeStepLength sub-stepping of numerical limit
+        /** atomicPhysics factor between actual timeStepLength in atomicPhysics sub-stepping of numerical limit
+         *
+         * @attention must be <= 1, otherwise solver is numerical unstable
+         */
         static constexpr float_X timeStepAlpha = 0.3_X;
 
-        // which probability approximation to use for the acceptance step
-        using ProbabilityApproximationFunctor
-            = picongpu::particles::atomicPhysics::ExponentialApproximationProbability;
+        /** maximum number of atomicPhysics sub-steps per PIC time step
+         *
+         * all rates above timeStepAlpha / maximumNumberSubStepsPerPICTimeStep will be solved in
+         *  rate fix point(equilibirium) approximation only.
+         *
+         * @attention This limit is currently only enforced for field ionization transitions! If the number of
+         * sub-steps required by collisional processes is higher than the limit set, AtomicPhysics will perform more
+         * subs-steps than set here.
+         */
+        static constexpr uint32_t maximumNumberSubStepsPerPICTimeStep = 200u;
     };
 
     /** atomicConfigNumber definition for argon ions

include/picongpu/param/atomicPhysics_Debug.param

@@ -140,8 +140,16 @@ namespace picongpu::atomicPhysics::debug
         constexpr bool CHECK_FOR_INVALID_TRANSITION_TYPE = false;
         //! @attention performance relevant
         constexpr bool CHECK_FOR_OVERFLOWS_IN_ACCUMULATON = false;
+        //! @attention performance relevant
+        constexpr bool CHECK_FOR_INFINITE_FIELD_IONIZATION_RATES = false;
     } // namespace kernel::chooseTransition
 
+    namespace kernel::chooseInstantTransition
+    {
+        //! @attention performance relevant
+        constexpr bool CHECK_FOR_INFINITE_FIELD_IONIZATION_RATES = false;
+    } // namespace kernel::chooseInstantTransition
+
     namespace kernel::recordSuggestedChanges
     {
         //! @attention only useful if using serial backend, no output unless compiling for cpu backend

include/picongpu/particles/atomicPhysics/InstantTransitionRateLimit.hpp

@@ -0,0 +1,45 @@
+/* Copyright 2024 Brian Marre
+ *
+ * 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
+
+// need unit system and normalization
+#include "picongpu/defines.hpp"
+#include "picongpu/particles/atomicPhysics/param.hpp"
+
+namespace picongpu::particles::atomicPhysics
+{
+    struct InstantTransitionRateLimit
+    {
+        /** get maximum of total state loss rate for inclusion in the time dependent rate equation solver
+         *
+         * @tparam T_ReturnType type and precision to use in the result
+         */
+        template<typename T_ReturnType>
+        static constexpr T_ReturnType get()
+        {
+            using picongpu::atomicPhysics::RateSolverParam;
+
+            // unit: unitless * unitless / unit_time = 1/unit_time
+            return static_cast<T_ReturnType>(
+                       RateSolverParam::timeStepAlpha * float_X(RateSolverParam::maximumNumberSubStepsPerPICTimeStep))
+                / picongpu::sim.pic.getDt<T_ReturnType>();
+        }
+    };
+} // namespace picongpu::particles::atomicPhysics

include/picongpu/particles/atomicPhysics/kernel/CalculateStepLength.kernel

@@ -91,10 +91,10 @@ namespace picongpu::particles::atomicPhysics::kernel
 
                     // 1/sim.unit.time()
                     float_X lossRateState = 0._X;
-                    if(rateCacheLossRate == 0._X)
+                    if((rateCacheLossRate <= 0._X))
                     {
-                        // no loss state, for example completely ionized state, set to neutral element
-                        // 1/sim.unit.time()
+                        /* no loss state, for example completely ionized state, set to neutral element
+                         * unit: 1/unit_time() */
                         lossRateState = picongpu::atomicPhysics::RateSolverParam::timeStepAlpha / sim.pic.getDt();
                     }
                     else
@@ -112,13 +112,16 @@ namespace picongpu::particles::atomicPhysics::kernel
                             lossRateState = picongpu::atomicPhysics::RateSolverParam::timeStepAlpha / sim.pic.getDt();
                         }
 
+                    // all state with field ionization loss rates above the limit have already been processed
+
                     // 1/(1/sim.unit.time()) = sim.unit.time()
                     float_X const speciesStateLimitTimeStep = 1._X / (lossRateState);
 
-                    // sim.unit.time()
+                    // unit: unit_time
                     alpaka::atomicMin(
                         worker.getAcc(),
-                        &timeStep, // sim.unit.time()
+                        // unit: unit_time
+                        &timeStep,
                         speciesStateLimitTimeStep * picongpu::atomicPhysics::RateSolverParam::timeStepAlpha,
                         ::alpaka::hierarchy::Threads{});
                 });

include/picongpu/particles/atomicPhysics/kernel/CheckForFieldEnergyOverSubscription.kernel

@@ -0,0 +1,122 @@
+/* Copyright 2024 Brian Marre
+ *
+ * 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/atomicPhysics/ConvertEnum.hpp"
+#include "picongpu/particles/atomicPhysics/KernelIndexation.hpp"
+#include "picongpu/particles/atomicPhysics/kernel/CalculateRejectionProbability.hpp"
+#include "picongpu/particles/atomicPhysics/localHelperFields/RejectionProbabilityCache_Cell.hpp"
+
+#include <pmacc/lockstep/ForEach.hpp>
+
+#include <cstdint>
+
+namespace picongpu::particles::atomicPhysics::kernel
+{
+    /** check all cells for overSubscription
+     *
+     * checks each superCell cell for fieldEnergyUse > fieldEnergy,
+     *  if yes mark cell as oversubscribed and stores rejection probability in rejection probability cache
+     *
+     * @tparam T_numberAtomicPhysicsIonSpecies specialization template parameter used to prevent compilation of all
+     *  atomicPhysics kernels if no atomic physics species is present.
+     */
+    template<uint32_t T_numberAtomicPhysicsIonSpecies>
+    struct CheckForFieldEnergyOverSubscriptionKernel
+    {
+        using VectorIdx = DataSpace<picongpu::simDim>;
+        /** call operator
+         *
+         * called by CheckForOverSubscription atomic physics sub-stage
+         *
+         * @param worker object containing the device and block information, passed by PMACC_KERNEL call
+         * @param areMapping mapping of blockIndex to block superCell index
+         * @param rngFactory factory for uniformly distributed random number generator
+         * @param eFieldBox deviceDataBox giving access to eField values for all local superCells
+         * @param fieldEnergyUseCacheBox deviceDataBox giving access to the field energy use cache for each local
+         * superCell
+         * @param sharedResourcesOverSubscribedBox deviceDataBox giving access to the local shared resources over
+         *  subscription switch for each local superCell
+         * @param rejectionProbabilityCacheCellBox deviceDataBox giving access to localRejectionProbabilityCache for
+         * all cells of all local superCells
+         */
+        template<
+            typename T_Worker,
+            typename T_AreaMapping,
+            typename T_TimeRemainingBox,
+            typename T_EFieldBox,
+            typename T_FieldEnergyUseCacheBox,
+            typename T_SharedRessourcesOverSubscribedBox,
+            typename T_RejectionProbabilityCacheCellDataBox>
+        HDINLINE void operator()(
+            T_Worker const& worker,
+            T_AreaMapping const areaMapping,
+            T_TimeRemainingBox const timeRemainingBox,
+            T_EFieldBox const eFieldBox,
+            T_FieldEnergyUseCacheBox const fieldEnergyUseCacheBox,
+            T_SharedRessourcesOverSubscribedBox sharedResourcesOverSubscribedBox,
+            T_RejectionProbabilityCacheCellDataBox rejectionProbabilityCacheCellBox) const
+        {
+            auto const superCellIdx = KernelIndexation::getSuperCellIndex(worker, areaMapping);
+            auto const superCellFieldIdx
+                = KernelIndexation::getSuperCellFieldIndexFromSuperCellIndex(areaMapping, superCellIdx);
+
+            auto const timeRemaining = timeRemainingBox(superCellFieldIdx);
+            if(timeRemaining <= 0._X)
+                return;
+
+            using FieldEnergyCache = typename T_FieldEnergyUseCacheBox::ValueType;
+            particles::atomicPhysics::localHelperFields::RejectionProbabilityCache_Cell<FieldEnergyCache::numberCells>&
+                rejectionProbabilityCacheCell
+                = rejectionProbabilityCacheCellBox(superCellFieldIdx);
+
+            FieldEnergyCache const& eFieldEnergyUseCache = fieldEnergyUseCacheBox(superCellFieldIdx);
+            VectorIdx const superCellCellOffset = superCellIdx * picongpu::SuperCellSize::toRT();
+            bool sharedResourcesOverSubscribed = false;
+
+            auto forEachCell = pmacc::lockstep::makeForEach<FieldEnergyCache::numberCells, T_Worker>(worker);
+            forEachCell(
+                [&worker,
+                 &superCellCellOffset,
+                 &eFieldBox,
+                 &eFieldEnergyUseCache,
+                 &rejectionProbabilityCacheCell,
+                 &sharedResourcesOverSubscribed](uint32_t const linearCellIndex)
+                {
+                    if(CalculateRejectionProbability::ofCell(
+                           linearCellIndex,
+                           superCellCellOffset,
+                           eFieldBox,
+                           eFieldEnergyUseCache,
+                           rejectionProbabilityCacheCell))
+                        sharedResourcesOverSubscribed = true;
+                });
+
+            uint32_t& flagField = sharedResourcesOverSubscribedBox(superCellFieldIdx);
+            // write out flag setting to device memory
+            alpaka::atomicOr(
+                worker.getAcc(),
+                &flagField,
+                u32(sharedResourcesOverSubscribed),
+                ::alpaka::hierarchy::Threads{});
+        }
+    };
+} // namespace picongpu::particles::atomicPhysics::kernel