Instant transition for states with high field ionization rate.

include/picongpu/param/atomicPhysics.param

@@ -55,6 +55,9 @@ namespace picongpu::atomicPhysics
         // 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 for fieldIonization
+        static constexpr uint32_t maximumNumberSubStepsFieldIonization = 2000;
     };
 
     /** atomicConfigNumber definition for argon ions

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

@@ -27,12 +27,16 @@
 
 #include <pmacc/lockstep/ForEach.hpp>
 
+#include <cstdint>
+
 namespace picongpu::particles::atomicPhysics::kernel
 {
     /** find atomicPhysics time step length kernel
      *
      * will find minimum time step length from the minimum stepLength of all atomic states
      *  of all species
+     *
+     * @tparam T_numberAtomicStates number of atomic states
      */
     template<uint32_t T_numberAtomicStates>
     struct CalculateStepLengthKernel

include/picongpu/particles/atomicPhysics/kernel/ChooseTransition_FieldBoundFree.kernel

@@ -176,15 +176,15 @@ namespace picongpu::particles::atomicPhysics::kernel
             DataBox<SharedBox<typename T_EField::ValueType, typename BlockArea::FullSuperCellSize, 0u>> eFieldCache
                 = CachedBox::create<0u, typename T_EField::ValueType>(worker, BlockArea());
 
-            // init E-Field cache
+            //      init
             auto const superCellSize = picongpu::SuperCellSize::toRT();
             DataSpace<picongpu::simDim> const superCellCellOffset = superCellIdx * superCellSize;
             auto fieldEBlockToCache = eFieldBox.shift(superCellCellOffset);
             pmacc::math::operation::Assign assign;
             auto collective = makeThreadCollective<BlockArea>();
             collective(worker, assign, eFieldCache, fieldEBlockToCache);
 
-            // wait for init to finish
+            //      wait for init to finish
             worker.sync();
 
             float_X const ionizationPotentialDepression

include/picongpu/particles/atomicPhysics/kernel/EFieldNormExtrema.hpp

@@ -0,0 +1,81 @@
+/* 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 simulation.param for normalisation and units, memory.param for SuperCellSize and dim.param for simDim
+#include "picongpu/defines.hpp"
+
+#include <pmacc/lockstep/ForEach.hpp>
+
+#include <limits>
+
+namespace picongpu::particles::atomicPhysics::kernel
+{
+    struct EFieldNormExtrema
+    {
+        //! find minimum and maximum value of the e E-Field Norm of the given superCell
+        template<typename T_Worker, typename T_EFieldDataBox, typename T_Type>
+        HDINLINE static void find(
+            T_Worker const& worker,
+            pmacc::DataSpace<picongpu::simDim> const superCellIdx,
+            T_EFieldDataBox const eFieldBox,
+            T_Type& minEFieldSuperCell,
+            T_Type& maxEFieldSuperCell)
+        {
+            auto onlyMaster = lockstep::makeMaster(worker);
+            onlyMaster(
+                [&maxEFieldSuperCell, &minEFieldSuperCell]()
+                {
+                    maxEFieldSuperCell = 0._X;
+                    minEFieldSuperCell = std::numeric_limits<float_X>::max();
+                });
+            worker.sync();
+
+            constexpr auto numberCellsInSuperCell
+                = pmacc::math::CT::volume<typename picongpu::SuperCellSize>::type::value;
+            DataSpace<picongpu::simDim> const superCellCellOffset = superCellIdx * picongpu::SuperCellSize::toRT();
+            auto forEachCell = pmacc::lockstep::makeForEach<numberCellsInSuperCell, T_Worker>(worker);
+
+            /// @todo switch to shared memory reduce, Brian Marre, 2024
+            forEachCell(
+                [&worker, &superCellCellOffset, &maxEFieldSuperCell, &minEFieldSuperCell, &eFieldBox](
+                    uint32_t const linearIdx)
+                {
+                    DataSpace<picongpu::simDim> const localCellIndex
+                        = pmacc::math::mapToND(picongpu::SuperCellSize::toRT(), static_cast<int>(linearIdx));
+                    DataSpace<picongpu::simDim> const cellIndex = localCellIndex + superCellCellOffset;
+
+                    auto const eFieldNorm = pmacc::math::l2norm(eFieldBox(cellIndex));
+
+                    alpaka::atomicMax(
+                        worker.getAcc(),
+                        // sim.unit.eField()
+                        &maxEFieldSuperCell,
+                        eFieldNorm);
+
+                    alpaka::atomicMin(
+                        worker.getAcc(),
+                        // sim.unit.eField()
+                        &minEFieldSuperCell,
+                        eFieldNorm);
+                });
+        }
+    };
+} // namespace picongpu::particles::atomicPhysics::kernel

include/picongpu/particles/atomicPhysics/kernel/FillLocalRateCache_BoundFree.kernel

@@ -26,6 +26,7 @@
 #include "picongpu/particles/atomicPhysics/enums/ADKLaserPolarization.hpp"
 #include "picongpu/particles/atomicPhysics/enums/ProcessClassGroup.hpp"
 #include "picongpu/particles/atomicPhysics/enums/TransitionOrdering.hpp"
+#include "picongpu/particles/atomicPhysics/kernel/EFieldNormExtrema.hpp"
 #include "picongpu/particles/atomicPhysics/rateCalculation/BoundFreeCollisionalTransitionRates.hpp"
 #include "picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp"
 
@@ -35,6 +36,7 @@
 #include <pmacc/static_assert.hpp>
 
 #include <cstdint>
+#include <limits>
 
 namespace picongpu::particles::atomicPhysics::kernel
 {
@@ -210,43 +212,20 @@ namespace picongpu::particles::atomicPhysics::kernel
             T_BoundFreeTransitionDataBox const boundFreeTransitionDataBox,
             float_X const ionizationPotentialDepression)
         {
-            // internal units
-            PMACC_SMEM(worker, maximumNormEFieldValueSuperCell, typename T_EFieldDataBox::ValueType::type);
+            // sim.unit.eField()
+            PMACC_SMEM(worker, minEFieldSuperCell, typename T_EFieldDataBox::ValueType::type);
+            // sim.unit.eField()
+            PMACC_SMEM(worker, maxEFieldSuperCell, typename T_EFieldDataBox::ValueType::type);
 
-            // find maximum field strength of superCell
-
-            /// @todo switch to shared memory reduce, Brian Marre, 2024
-            constexpr auto numberCellsInSuperCell
-                = pmacc::math::CT::volume<typename picongpu::SuperCellSize>::type::value;
-            auto forEachCell = pmacc::lockstep::makeForEach<numberCellsInSuperCell, T_Worker>(worker);
-
-            DataSpace<picongpu::simDim> const superCellCellOffset = superCellIdx * picongpu::SuperCellSize::toRT();
-
-            /// @todo switch to shared memory reduce, Brian Marre, 2024
-            forEachCell(
-                [&worker, &superCellCellOffset, &maximumNormEFieldValueSuperCell, &eFieldBox](uint32_t const linearIdx)
-                {
-                    DataSpace<picongpu::simDim> const localCellIndex
-                        = pmacc::math::mapToND(picongpu::SuperCellSize::toRT(), static_cast<int>(linearIdx));
-                    DataSpace<picongpu::simDim> const cellIndex = localCellIndex + superCellCellOffset;
-
-                    auto const eFieldNorm = pmacc::math::l2norm(eFieldBox(cellIndex));
-
-                    alpaka::atomicMax(
-                        worker.getAcc(),
-                        // internal units
-                        &maximumNormEFieldValueSuperCell,
-                        eFieldNorm);
-                });
+            EFieldNormExtrema::find(worker, superCellIdx, eFieldBox, minEFieldSuperCell, maxEFieldSuperCell);
             worker.sync();
 
             // calculate ADK field ionization rate, for each atomic state and maximum field value
             auto forEachAtomicState = pmacc::lockstep::makeForEach<T_numberAtomicStates, T_Worker>(worker);
-
-            //      accumulate rates of possible transitions by state
             forEachAtomicState(
                 [&ionizationPotentialDepression,
-                 &maximumNormEFieldValueSuperCell,
+                 &maxEFieldSuperCell,
+                 &minEFieldSuperCell,
                  &rateCache,
                  &numberTransitionsDataBox,
                  &startIndexDataBox,
@@ -269,15 +248,18 @@ namespace picongpu::particles::atomicPhysics::kernel
                         uint32_t const transitionCollectionIndex = offset + transitionID;
 
                         // 1/picongpu::sim.unit.time()
-                        sumRateTransitions += atomicPhysics::rateCalculation::
-                            BoundFreeFieldTransitionRates<T_ADKLaserPolarization>::template rateADKFieldIonization(
-                                maximumNormEFieldValueSuperCell,
-                                ionizationPotentialDepression,
-                                transitionCollectionIndex,
-                                chargeStateDataDataBox,
-                                atomicStateDataDataBox,
-                                boundFreeTransitionDataBox);
+                        sumRateTransitions
+                            += atomicPhysics::rateCalculation::BoundFreeFieldTransitionRates<T_ADKLaserPolarization>::
+                                template maximumRateADKFieldIonization(
+                                    minEFieldSuperCell,
+                                    maxEFieldSuperCell,
+                                    ionizationPotentialDepression,
+                                    transitionCollectionIndex,
+                                    chargeStateDataDataBox,
+                                    atomicStateDataDataBox,
+                                    boundFreeTransitionDataBox);
                     }
+
                     rateCache.template add<s_enums::ChooseTransitionGroup::fieldBoundFreeUpward>(
                         atomicStateCollectionIndex,
                         sumRateTransitions);

include/picongpu/particles/atomicPhysics/kernel/RecordChanges.kernel

@@ -125,10 +125,14 @@ namespace picongpu::particles::atomicPhysics::kernel
                 return;
 
             // get histogram for current superCell
-            T_Histogram& electronHistogram = localElectronHistogramDataBox(superCellFieldIdx);
+            [[maybe_unused]] T_Histogram& electronHistogram = localElectronHistogramDataBox(superCellFieldIdx);
 
             // eV
             [[maybe_unused]] float_X ionizationPotentialDepression = 0._X;
+
+            /* field ionization is not energy conserving currently
+             *-> no need to calculate energy used -> no need to calculate IPD */
+            /// @todo implement energy conservation in field ionization, Brian Marre, 2024
             if constexpr(isIonizing && isCollisional)
                 ionizationPotentialDepression
                     = ionizationPotentialDepression::PassIPDInputs::template calculateIPD<T_IPDModel>(