Merge pull request #5245 from BrianMarre/topic-bugFixIPDIonization

bug fix ipd ionization
ComputationalRadiationPhysics · Jan 10, 2025 · 74990f4 · 74990f4
2 parents 648c864 + 6eb0d91
commit 74990f4
Show file tree

Hide file tree

Showing 7 changed files with 31 additions and 17 deletions.
diff --git a/include/picongpu/particles/atomicPhysics/DeltaEnergyTransition.hpp b/include/picongpu/particles/atomicPhysics/DeltaEnergyTransition.hpp
@@ -122,8 +122,8 @@ namespace picongpu::particles::atomicPhysics
 
         /** get deltaEnergy of transition
          *
-         * DeltaEnergy is defined as energy(UpperState) - energy(lowerState) [+ ionizationEnergy],
-         *  with lower and upper state as given in charge state box
+         * DeltaEnergy is defined as energy(upperState) - energy(lowerState) [+ ionizationEnergy],
+         *  with lower and upper state specified by the given transitionCollectionIndex
          *
          * @param atomicStateBox deviceDataBox giving access to atomic state property data
          * @param transitionBox deviceDataBox giving access to transition property data,

diff --git a/include/picongpu/particles/atomicPhysics/atomicData/AtomicStateData.hpp b/include/picongpu/particles/atomicPhysics/atomicData/AtomicStateData.hpp
@@ -40,7 +40,7 @@ namespace picongpu::particles::atomicPhysics::atomicData
      *
      * @tparam T_Number dataType used for number storage, typically uint32_t
      * @tparam T_Value dataType used for value storage, typically float_X
-     * @tparam T_ConfigNumber dataType used for storage of configNumber of atomic states
+     * @tparam T_ConfigNumber dataType used for conversion from/to configNumber to/from super configuration
      * @tparam T_Multiplicity dataType used for T_Multiplicity storage, typically uint64_t
      *
      * @attention ConfigNumber specifies the number of a state as defined by the configNumber

diff --git a/include/picongpu/particles/atomicPhysics/debug/PrintAtomicPhysicsIonToConsole.hpp b/include/picongpu/particles/atomicPhysics/debug/PrintAtomicPhysicsIonToConsole.hpp
@@ -29,7 +29,7 @@ namespace picongpu::particles::atomicPhysics::debug
 {
     /** debug only, write atomicPhysics attributes to console
      *
-     * @attention only creates ouptut if atomicPhysics debug setting CPU_OUTPUT_ACTIVE == True
+     * @attention only creates output if compiling for debug backend
      * @attention only useful if compiling serial backend
      */
     struct PrintAtomicPhysicsIonToConsole

diff --git a/include/picongpu/particles/atomicPhysics/debug/TestRateCalculation.hpp b/include/picongpu/particles/atomicPhysics/debug/TestRateCalculation.hpp
@@ -406,18 +406,18 @@ namespace picongpu::particles::atomicPhysics::debug
         //! @return true =^= test passed
         bool testADKIonizationRate() const
         {
-            // 1/s
+            // unit: 1/s
             float_64 const correctRate = 6.391666527e+9 * 1 / 3.3e-17;
 
-            // sim.unit.eField()
+            // unit: unit_eField
             float_X const eFieldNorm = 0.0126 * sim.atomicUnit.eField() / sim.unit.eField();
 
-            // eV
+            // unit: eV
             float_X const ipd = 0._X;
 
             constexpr auto laserPolarization
                 = picongpu::particles::atomicPhysics::enums::ADKLaserPolarization::linearPolarization;
-            // 1/s
+            // unit: 1/s
             float_64 const rate
                 = static_cast<float_64>(
                       rateCalculation::BoundFreeFieldTransitionRates<laserPolarization>::rateADKFieldIonization<

diff --git a/...pu/particles/atomicPhysics/ionizationPotentialDepression/kernel/ApplyIPDIonization.kernel b/...pu/particles/atomicPhysics/ionizationPotentialDepression/kernel/ApplyIPDIonization.kernel
@@ -201,12 +201,23 @@ namespace picongpu::particles::atomicPhysics::ionizationPotentialDepression::ker
                 ipd += BarrierSupressionIonization::getIPD(screenedCharge, eFieldNormAU);
             }
 
-            // eV
-            float_X const ipdIonizationEnergy = ionizationEnergy - ipd;
+            bool const stateIsUnbound = ((ionizationEnergy - ipd) <= 0._X);
+            bool const stateIsGenerallyUnbound = (ionizationEnergy <= 0._X);
 
-            if(ipdIonizationEnergy < 0._X)
+            /** @details states that are unbound without an IPD contribution must relax via auto-ionization, electronic
+             *  collisional ionization or deexcitation channels in the regular rate solver since IPD lacks the energy
+             *  for ionization.
+             *
+             *  These states are typically low shell hole states, IPD-Ionization will not relax these states until
+             *  very high charge states, causing:
+             *  - numerical energy creation as the IPD actually lacks the energy required for ionization
+             *  - instant ionization cascades, which are not consistent with the IPD equilibrium description
+             *
+             * @todo improve by developing a non equlibirium IPD description, Brian Marre, 2025
+             */
+            if(!stateIsGenerallyUnbound && stateIsUnbound)
             {
-                /* we only update the atomic state since IPDIonization is not a regular transition and does not use
+                /* we only update the atomic state since IPD-Ionization is not a regular transition and does not use
                  *   shared resources */
 
                 // update ion atomic state

diff --git a/include/picongpu/particles/atomicPhysics/kernel/FillRateCache_Autonomous.kernel b/include/picongpu/particles/atomicPhysics/kernel/FillRateCache_Autonomous.kernel
@@ -133,13 +133,13 @@ namespace picongpu::particles::atomicPhysics::kernel
 
                         for(uint32_t i = u32(0u); i < numberTransitionsDown; i++)
                         {
+                            uint32_t const transitionCollectionIndex = offset + i;
                             // 1/sim.unit.time()
                             rateCache.template add<atomicPhysics::enums::ChooseTransitionGroup::autonomousDownward>(
                                 atomicStateCollectionIndex,
                                 picongpu::particles::atomicPhysics::rateCalculation::AutonomousTransitionRates::
                                     template rateAutonomousIonization<T_AutonomousTransitionDataBox>(
-                                        // transitionCollectionIndex
-                                        offset + i,
+                                        transitionCollectionIndex,
                                         autonomousTransitionDataBox));
                         }
                     });

diff --git a/.../picongpu/particles/atomicPhysics/rateCalculation/BoundFreeCollisionalTransitionRates.hpp b/.../picongpu/particles/atomicPhysics/rateCalculation/BoundFreeCollisionalTransitionRates.hpp
@@ -169,9 +169,12 @@ namespace picongpu::particles::atomicPhysics::rateCalculation
                 /// @todo replace with screenedCharge(upperStateChargeState)?, Brian Marre, 2023
                 float_X const screenedCharge = chargeStateDataBox.screenedCharge(lowerStateChargeState) - 1._X; // [e]
 
-                float_X const U = energyElectron / energyDifference; // unitless
-                float_X const beta = betaFactor(screenedCharge); // unitless
-                float_64 const w = wFactor(U, beta); // unitless
+                // unitless
+                float_X const U = energyElectron / energyDifference;
+                // unitless
+                float_X const beta = betaFactor(screenedCharge);
+                // unitless
+                float_64 const w = wFactor(U, beta);
                 float_X const crossSection = static_cast<float_X>(
                     scalingConstant * combinatorialFactor
                     / static_cast<float_64>(pmacc::math::cPow(energyDifference, 2u)) / static_cast<float_64>(U)