diff --git a/include/picongpu/param/atomicPhysics.param b/include/picongpu/param/atomicPhysics.param index d2478cc4b8..9b101afadc 100644 --- a/include/picongpu/param/atomicPhysics.param +++ b/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 diff --git a/include/picongpu/param/atomicPhysics_Debug.param b/include/picongpu/param/atomicPhysics_Debug.param index baaf8cf8a1..4d556dbb5f 100644 --- a/include/picongpu/param/atomicPhysics_Debug.param +++ b/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 diff --git a/include/picongpu/particles/atomicPhysics/InstantTransitionRateLimit.hpp b/include/picongpu/particles/atomicPhysics/InstantTransitionRateLimit.hpp new file mode 100644 index 0000000000..323fdbfe0d --- /dev/null +++ b/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 . + */ + +#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 + static constexpr T_ReturnType get() + { + using picongpu::atomicPhysics::RateSolverParam; + + // unit: unitless * unitless / unit_time = 1/unit_time + return static_cast( + RateSolverParam::timeStepAlpha * float_X(RateSolverParam::maximumNumberSubStepsPerPICTimeStep)) + / picongpu::sim.pic.getDt(); + } + }; +} // namespace picongpu::particles::atomicPhysics diff --git a/include/picongpu/particles/atomicPhysics/kernel/CalculateStepLength.kernel b/include/picongpu/particles/atomicPhysics/kernel/CalculateStepLength.kernel index b594fa2e8d..417992c45d 100644 --- a/include/picongpu/particles/atomicPhysics/kernel/CalculateStepLength.kernel +++ b/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{}); }); diff --git a/include/picongpu/particles/atomicPhysics/kernel/CheckForFieldEnergyOverSubscription.kernel b/include/picongpu/particles/atomicPhysics/kernel/CheckForFieldEnergyOverSubscription.kernel new file mode 100644 index 0000000000..a53c4d13cf --- /dev/null +++ b/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 . + */ + +#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 + +#include + +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 + struct CheckForFieldEnergyOverSubscriptionKernel + { + using VectorIdx = DataSpace; + /** 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& + rejectionProbabilityCacheCell + = rejectionProbabilityCacheCellBox(superCellFieldIdx); + + FieldEnergyCache const& eFieldEnergyUseCache = fieldEnergyUseCacheBox(superCellFieldIdx); + VectorIdx const superCellCellOffset = superCellIdx * picongpu::SuperCellSize::toRT(); + bool sharedResourcesOverSubscribed = false; + + auto forEachCell = pmacc::lockstep::makeForEach(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 diff --git a/include/picongpu/particles/atomicPhysics/kernel/ChooseInstantNonReversibleTransition.kernel b/include/picongpu/particles/atomicPhysics/kernel/ChooseInstantNonReversibleTransition.kernel new file mode 100644 index 0000000000..086d9aefae --- /dev/null +++ b/include/picongpu/particles/atomicPhysics/kernel/ChooseInstantNonReversibleTransition.kernel @@ -0,0 +1,311 @@ +/* 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 . + */ + +#pragma once + +// need simulation.param for normalisation and units, memory.param for SuperCellSize and dim.param for simDim +#include "picongpu/defines.hpp" +#include "picongpu/particles/atomicPhysics/CheckSetOfAtomicDataBoxes.hpp" +#include "picongpu/particles/atomicPhysics/ConvertEnum.hpp" +#include "picongpu/particles/atomicPhysics/EFieldCache.hpp" +#include "picongpu/particles/atomicPhysics/InstantTransitionRateLimit.hpp" +#include "picongpu/particles/atomicPhysics/KernelIndexation.hpp" +#include "picongpu/particles/atomicPhysics/MinimumAndMaximumEFieldNormOfSuperCell.hpp" +#include "picongpu/particles/atomicPhysics/enums/ADKLaserPolarization.hpp" +#include "picongpu/particles/atomicPhysics/enums/ProcessClass.hpp" +#include "picongpu/particles/atomicPhysics/enums/TransitionDirection.hpp" +#include "picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp" + +#include +#include +#include + +#include + +namespace picongpu::particles::atomicPhysics::kernel +{ + namespace s_enums = picongpu::particles::atomicPhysics::enums; + + template + struct ChooseInstantNonReversibleTransitionKernel + { + using VectorIdx = pmacc::DataSpace; + + /** call operator + * + * called by ChooseInstantNonReversibleTransition atomic physics sub-stage + * + * @param worker object containing the device and block information, passed by PMACC_KERNEL call + * @param areaMapping mapping of blockIndex to block superCell index + * @param rngFactoryFloat factory for uniformly distributed random number generator, for float_X [0,1) + * @param chargeStateBox deviceDataBox giving access to charge state property data + * @param atomicStateBox deviceDataBox giving access to atomic state property data + * @param transitionBox deviceDataBox giving access to bound free transition data + * @param numberTransitionsBox deviceDataBox giving access to the number of + * bound-free transitions for each atomic state + * @param startIndexBox deviceDataBox giving access to the start index of each + * @param timeRemainingBox deviceDataBox giving access to the local time remaining of all local superCells + * @param eFieldBox deviceDataBox giving access to the device local E-Field values + * @param ionBox deviceDataBox giving access to the ion frames of all device local superCells + * @param ipdInput everything required by T_IPDModel to calculate the IonizationPotentialDepression, + * passed by T_IPDModel::callKernelWithIPDInput + * + * @attention assumes that the accepted ion attribute of the species has been reset before the kernel call + * @attention assumes that the foundUnbound-field has been reset before the kernel call + */ + template< + typename T_Worker, + typename T_AreaMapping, + typename T_RngGeneratorFactoryFloat, + typename T_LocalTimeRemainingBox, + typename T_FoundUnboundBox, + typename T_ChargeStateDataDataBox, + typename T_AtomicStateDataDataBox, + typename T_AtomicStateBoundFreeStartIndexBlockDataBox, + typename T_AtomicStateBoundFreeNumberTransitionsDataBox, + typename T_BoundFreeTransitionDataBox, + typename T_EFieldDataBox, + typename T_IonBox, + typename... T_IPDInput> + HDINLINE void operator()( + T_Worker const& worker, + T_AreaMapping const areaMapping, + T_RngGeneratorFactoryFloat rngFactoryFloat, + T_ChargeStateDataDataBox const chargeStateBox, + T_AtomicStateDataDataBox const atomicStateBox, + T_AtomicStateBoundFreeStartIndexBlockDataBox const startIndexBox, + T_AtomicStateBoundFreeNumberTransitionsDataBox const numberTransitionsBox, + T_BoundFreeTransitionDataBox const transitionBox, + T_LocalTimeRemainingBox const timeRemainingBox, + T_FoundUnboundBox foundUnboundBox, + T_EFieldDataBox const eFieldBox, + T_IonBox ionBox, + T_IPDInput const... ipdInput) const + { + PMACC_CASSERT(CheckSetOfAtomicDataBoxes::areBoundFreeAndDirection< + s_enums::TransitionDirection::upward, + T_AtomicStateBoundFreeStartIndexBlockDataBox, + T_AtomicStateBoundFreeNumberTransitionsDataBox, + T_BoundFreeTransitionDataBox>()); + + VectorIdx const superCellIdx = KernelIndexation::getSuperCellIndex(worker, areaMapping); + VectorIdx const superCellFieldIdx + = KernelIndexation::getSuperCellFieldIndexFromSuperCellIndex(areaMapping, superCellIdx); + + auto const timeRemaining = timeRemainingBox(superCellFieldIdx); + auto forEachLocalIonBoxEntry = pmacc::particles::algorithm::acc::makeForEach(worker, ionBox, superCellIdx); + + if((timeRemaining <= 0._X) || (!forEachLocalIonBoxEntry.hasParticles())) + return; + + // create cache for which states to apply kernel to + PMACC_SMEM(worker, stateHasInstantNonReversibleTransition, memory::Array); + // unit: unit_eField + PMACC_SMEM(worker, minEFieldNormSuperCell, typename T_EFieldDataBox::ValueType::type); + // unit: unit_eField + PMACC_SMEM(worker, maxEFieldNormSuperCell, typename T_EFieldDataBox::ValueType::type); + + // init shared memory + //!@{ + /// @todo write variant using the eField cache to find minimum and maximum using, Brian Marre, 2024 + MinimumAndMaximumEFieldNormOfSuperCell::find( + worker, + superCellIdx, + eFieldBox, + minEFieldNormSuperCell, + maxEFieldNormSuperCell); + auto forEachAtomicState = pmacc::lockstep::makeForEach(worker); + forEachAtomicState([&stateHasInstantNonReversibleTransition](uint32_t const stateCollectionIndex) + { stateHasInstantNonReversibleTransition[stateCollectionIndex] = false; }); + worker.sync(); + //!@} + + // check presence of atomic states + forEachLocalIonBoxEntry( + [&stateHasInstantNonReversibleTransition](T_Worker const& worker, auto& ion) + { + uint32_t const stateCollectionIndex = ion[atomicStateCollectionIndex_]; + stateHasInstantNonReversibleTransition[stateCollectionIndex] = true; + }); + worker.sync(); + + float_X const ionizationPotentialDepression + = T_IPDModel::template calculateIPD( + superCellFieldIdx, + ipdInput...); + + // mask states with maximum total loss rate below rate limit + forEachAtomicState( + [&ionizationPotentialDepression, + &minEFieldNormSuperCell, + &maxEFieldNormSuperCell, + &stateHasInstantNonReversibleTransition, + &numberTransitionsBox, + &startIndexBox, + &chargeStateBox, + &atomicStateBox, + &transitionBox](uint32_t const stateCollectionIndex) + { + // check if atomic state present at all + if(!stateHasInstantNonReversibleTransition[stateCollectionIndex]) + return; + + uint32_t const numberTransitionsUp + = numberTransitionsBox.numberOfTransitionsUp(stateCollectionIndex); + uint32_t const offset = startIndexBox.startIndexBlockTransitionsUp(stateCollectionIndex); + + // unit: 1/unit_time + float_64 sumRateTransitions = 0.0; + for(uint32_t transitionID = u32(0u); transitionID < numberTransitionsUp; ++transitionID) + { + uint32_t const transitionCollectionIndex = offset + transitionID; + + // unit: 1/unit_time + sumRateTransitions + += atomicPhysics::rateCalculation::BoundFreeFieldTransitionRates:: + template maximumRateADKFieldIonization( + minEFieldNormSuperCell, + maxEFieldNormSuperCell, + ionizationPotentialDepression, + transitionCollectionIndex, + chargeStateBox, + atomicStateBox, + transitionBox); + } + + // unit: 1/unit_time + constexpr auto rateLimit = InstantTransitionRateLimit::get(); + + if(sumRateTransitions <= rateLimit) + stateHasInstantNonReversibleTransition[stateCollectionIndex] = false; + }); + // no need to snyc threads since we will sync after filling eField cache anyway + + auto eFieldCache = EFieldCache::get<__COUNTER__>(worker, superCellIdx, eFieldBox); + worker.sync(); + + auto rngGeneratorFloat = rngFactoryFloat(worker, superCellFieldIdx); + + bool foundInstantTransitionIon = false; + // choose instant transition for each particle in state with total transition rate above limit + forEachLocalIonBoxEntry( + [&stateHasInstantNonReversibleTransition, + &foundInstantTransitionIon, + &ionizationPotentialDepression, + &eFieldCache, + &rngGeneratorFloat, + &numberTransitionsBox, + &startIndexBox, + &chargeStateBox, + &atomicStateBox, + &transitionBox](T_Worker const& worker, auto& ion) + { + auto const stateCollectionIndex = ion[atomicStateCollectionIndex_]; + auto const linearCellIndex = ion[localCellIdx_]; + + if(!stateHasInstantNonReversibleTransition[stateCollectionIndex]) + return; + + VectorIdx const localCellIndex + = pmacc::math::mapToND(picongpu::SuperCellSize::toRT(), static_cast(linearCellIndex)); + + // unit: unit_eField + float_X const eFieldNormCell = pmacc::math::l2norm(eFieldCache(localCellIndex)); + + uint32_t const numberTransitions + = numberTransitionsBox.numberOfTransitionsUp(stateCollectionIndex); + uint32_t const startIndexTransitionBlock + = startIndexBox.startIndexBlockTransitionsUp(stateCollectionIndex); + + float_X const r = rngGeneratorFloat(); + + float_64 particleTotalTransitionRate = 0.; + for(uint32_t transitionID = 0u; transitionID < numberTransitions; ++transitionID) + { + uint32_t const transitionCollectionIndex = transitionID + startIndexTransitionBlock; + + auto const transitionRate + = atomicPhysics::rateCalculation::BoundFreeFieldTransitionRates:: + template rateADKFieldIonization( + eFieldNormCell, + ionizationPotentialDepression, + transitionCollectionIndex, + chargeStateBox, + atomicStateBox, + transitionBox); + + if constexpr(picongpu::atomicPhysics::debug::kernel::chooseInstantTransition:: + CHECK_FOR_INFINITE_FIELD_IONIZATION_RATES) + if(transitionRate < 0._X) + printf("atomicPhysics ERROR: encountered infinite transition rate in " + "ChooseInstantTransition kernel"); + + // unit: 1/unit_time + particleTotalTransitionRate += transitionRate; + } + + constexpr auto rateLimit = InstantTransitionRateLimit::get(); + bool const particleTotalLossRateBelowRateLimit = (particleTotalTransitionRate <= rateLimit); + + if(particleTotalLossRateBelowRateLimit) + return; + + foundInstantTransitionIon = true; + + // choose transition for instant transition ion + float_X cumSum = 0._X; + for(uint32_t transitionID = 0u; transitionID < numberTransitions; ++transitionID) + { + uint32_t const transitionCollectionIndex = transitionID + startIndexTransitionBlock; + + // unit: 1/unit_time + cumSum += static_cast( + atomicPhysics::rateCalculation::BoundFreeFieldTransitionRates:: + template rateADKFieldIonization( + eFieldNormCell, + ionizationPotentialDepression, + transitionCollectionIndex, + chargeStateBox, + atomicStateBox, + transitionBox) + / particleTotalTransitionRate); + + // inclusive limit, to make sure that r==1 is assigned a transition + if(r <= cumSum) + { + // found chosen transition + ion[processClass_] = u8(s_enums::ProcessClass::fieldIonization); + ion[transitionIndex_] = transitionCollectionIndex; + /* field ionizations are not bin based therefore we do not set a bin, and old values are + * ignored */ + // we set the accepted flag to allow easy resource use accounting in a later kernel call + ion[accepted_] = true; + return; + } + } + }); + + alpaka::atomicOr( + worker.getAcc(), + &foundUnboundBox(superCellFieldIdx), + u32(foundInstantTransitionIon), + ::alpaka::hierarchy::Threads{}); + } + }; +} // namespace picongpu::particles::atomicPhysics::kernel diff --git a/include/picongpu/particles/atomicPhysics/kernel/ChooseTransitionGroup.kernel b/include/picongpu/particles/atomicPhysics/kernel/ChooseTransitionGroup.kernel index 0a1c4abb13..b4177c2a84 100644 --- a/include/picongpu/particles/atomicPhysics/kernel/ChooseTransitionGroup.kernel +++ b/include/picongpu/particles/atomicPhysics/kernel/ChooseTransitionGroup.kernel @@ -164,7 +164,6 @@ namespace picongpu::particles::atomicPhysics::kernel if(r < cumSum) { // only in between storage - // printf("ID: %lu, %i\n", ion[particleId_], chooseTransitionGroupID); ion[transitionIndex_] = chooseTransitionGroupID; break; } diff --git a/include/picongpu/particles/atomicPhysics/kernel/ChooseTransition_FieldBoundFree.kernel b/include/picongpu/particles/atomicPhysics/kernel/ChooseTransition_FieldBoundFree.kernel index 2e8b930ac1..e004b23e96 100644 --- a/include/picongpu/particles/atomicPhysics/kernel/ChooseTransition_FieldBoundFree.kernel +++ b/include/picongpu/particles/atomicPhysics/kernel/ChooseTransition_FieldBoundFree.kernel @@ -183,6 +183,12 @@ namespace picongpu::particles::atomicPhysics::kernel atomicStateDataDataBox, transitionDataBox); + if constexpr(picongpu::atomicPhysics::debug::kernel::chooseTransition:: + CHECK_FOR_INFINITE_FIELD_IONIZATION_RATES) + if(rateTransition < 0._X) + printf("atomicPhysics ERROR: encountered infinite transition rate in non instant " + "transition ion"); + cumSum += rateTransition / rateCache.rate( u32(s_enums::ChooseTransitionGroup::fieldBoundFreeUpward), diff --git a/include/picongpu/particles/atomicPhysics/kernel/FillRateCache_BoundFree.kernel b/include/picongpu/particles/atomicPhysics/kernel/FillRateCache_BoundFree.kernel index f558358ffa..87effe911b 100644 --- a/include/picongpu/particles/atomicPhysics/kernel/FillRateCache_BoundFree.kernel +++ b/include/picongpu/particles/atomicPhysics/kernel/FillRateCache_BoundFree.kernel @@ -23,6 +23,7 @@ #include "picongpu/defines.hpp" #include "picongpu/particles/atomicPhysics/CheckSetOfAtomicDataBoxes.hpp" #include "picongpu/particles/atomicPhysics/ConvertEnum.hpp" +#include "picongpu/particles/atomicPhysics/InstantTransitionRateLimit.hpp" #include "picongpu/particles/atomicPhysics/KernelIndexation.hpp" #include "picongpu/particles/atomicPhysics/MinimumAndMaximumEFieldNormOfSuperCell.hpp" #include "picongpu/particles/atomicPhysics/electronDistribution/CachedHistogram.hpp" @@ -231,9 +232,8 @@ namespace picongpu::particles::atomicPhysics::kernel { uint32_t const transitionCollectionIndex = offset + transitionID; - // unit: 1/unit_time - sumRateTransitions - += atomicPhysics::rateCalculation::BoundFreeFieldTransitionRates:: + auto const transitionRate + = atomicPhysics::rateCalculation::BoundFreeFieldTransitionRates:: template maximumRateADKFieldIonization( minEFieldNormSuperCell, maxEFieldNormSuperCell, @@ -242,6 +242,23 @@ namespace picongpu::particles::atomicPhysics::kernel chargeStateDataDataBox, atomicStateDataDataBox, boundFreeTransitionDataBox); + + if(transitionRate > 0._X) + { + // unit: 1/unit_time + sumRateTransitions += transitionRate; + } + } + + /* all particles with rates above the instant transition limit are processed in the + * chooseInstantRateTransitions kernel therefore the remaining particles always have a state loss + * rate limit below or equal to the limit*/ + + constexpr auto instantTransitionRateLimit = InstantTransitionRateLimit::get(); + bool const stateWithInstantTransitions = (sumRateTransitions > instantTransitionRateLimit); + if(stateWithInstantTransitions) + { + sumRateTransitions = instantTransitionRateLimit; } rateCache.template add( diff --git a/include/picongpu/particles/atomicPhysics/kernel/RecordSuggestedFieldEnergyUse.kernel b/include/picongpu/particles/atomicPhysics/kernel/RecordSuggestedFieldEnergyUse.kernel new file mode 100644 index 0000000000..37577b54dc --- /dev/null +++ b/include/picongpu/particles/atomicPhysics/kernel/RecordSuggestedFieldEnergyUse.kernel @@ -0,0 +1,120 @@ +/* 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/particles/atomicPhysics/DeltaEnergyTransition.hpp" +#include "picongpu/particles/atomicPhysics/KernelIndexation.hpp" +#include "picongpu/particles/atomicPhysics/enums/IsProcess.hpp" +#include "picongpu/particles/atomicPhysics/enums/ProcessClassGroup.hpp" + +#include + +#include + +namespace picongpu::particles::atomicPhysics::kernel +{ + /** record used electric field energy from a cell for each accepted transition + * + * @tparam T_IPDModel ionization potential depression model to use + */ + template + struct RecordSuggestedFieldEnergyUseKernel + { + /** call operator + * + * called by RecordSuggestedFieldEnergyUsed atomic physics sub-stage + * + * @param worker object containing the device and block information, passed by PMACC_KERNEL call + * @param areaMapping mapping of blockIndex to block superCell index + * @param chargeStateBox deviceDataBox giving access to charge state property data + * @param atomicStateBox deviceDataBox giving access to atomic state property data + * @param boundFreeTransitionBox deviceDataBox giving access to bound free transition data + * @param timeRemainingBox deviceDataBox giving access to the local time remaining of all local superCells + * @param fieldEnergyUseCacheBox deviceDataBox giving access to the local field energy use of all cells of all + * local superCells + * @param ionBox deviceDataBox giving access to the particle frames of all local superCells + * @param ipdInput everything required by T_IPDModel to calculate the IonizationPotentialDepression, + * passed by T_IPDModel::callKernelWithIPDInput + */ + template< + typename T_Worker, + typename T_AreaMapping, + typename T_ChargeStateDataDataBox, + typename T_AtomicStateDataDataBox, + typename T_BoundFreeTransitionDataBox, + typename T_LocalTimeRemainingBox, + typename T_FieldEnergyUseCacheBox, + typename T_IonBox, + typename... T_IPDInput> + HDINLINE void operator()( + T_Worker const& worker, + T_AreaMapping const areaMapping, + T_ChargeStateDataDataBox const chargeStateBox, + T_AtomicStateDataDataBox const atomicStateBox, + T_BoundFreeTransitionDataBox const boundFreeTransitionBox, + T_LocalTimeRemainingBox const timeRemainingBox, + T_FieldEnergyUseCacheBox fieldEnergyUseCacheBox, + T_IonBox const ionBox, + T_IPDInput const... ipdInput) const + { + auto const superCellIdx = KernelIndexation::getSuperCellIndex(worker, areaMapping); + auto const superCellFieldIdx + = KernelIndexation::getSuperCellFieldIndexFromSuperCellIndex(areaMapping, superCellIdx); + + auto const timeRemaining = timeRemainingBox(superCellFieldIdx); + auto forEachLocalIonBoxEntry = pmacc::particles::algorithm::acc::makeForEach(worker, ionBox, superCellIdx); + + // end kernel if no particles or superCell already finished + if((timeRemaining <= 0._X) || (!forEachLocalIonBoxEntry.hasParticles())) + return; + + auto& fieldEnergyUseCache = fieldEnergyUseCacheBox(superCellFieldIdx); + + float_X const ionizationPotentialDepression + = T_IPDModel::template calculateIPD( + superCellFieldIdx, + ipdInput...); + + forEachLocalIonBoxEntry( + [&](T_Worker const& worker, auto& ion) + { + if(!ion[accepted_]) + return; + + auto const weight = ion[weighting_]; + uint32_t const linearCellIdx = ion[localCellIdx_]; + uint32_t const transitionCollectionIndex = ion[transitionIndex_]; + + // eV, weighted + float_X const energyUsed = DeltaEnergyTransition::get( + transitionCollectionIndex, + atomicStateBox, + boundFreeTransitionBox, + ionizationPotentialDepression, + chargeStateBox) + * weight; + + // uses atomics + fieldEnergyUseCache.add(worker, linearCellIdx, energyUsed); + }); + } + }; +} // namespace picongpu::particles::atomicPhysics::kernel diff --git a/include/picongpu/particles/atomicPhysics/kernel/UpdateIonAtomicState.kernel b/include/picongpu/particles/atomicPhysics/kernel/UpdateIonAtomicState.kernel new file mode 100644 index 0000000000..8191f30cb9 --- /dev/null +++ b/include/picongpu/particles/atomicPhysics/kernel/UpdateIonAtomicState.kernel @@ -0,0 +1,116 @@ +/* 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 . + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/particles/atomicPhysics/CheckSetOfAtomicDataBoxes.hpp" +#include "picongpu/particles/atomicPhysics/ConvertEnum.hpp" +#include "picongpu/particles/atomicPhysics/KernelIndexation.hpp" +#include "picongpu/particles/atomicPhysics/SetAtomicState.hpp" +#include "picongpu/particles/atomicPhysics/enums/ProcessClass.hpp" +#include "picongpu/particles/atomicPhysics/enums/TransitionDirection.hpp" +#include "picongpu/particles/atomicPhysics/ionizationPotentialDepression/PassIPDInputs.hpp" + +#include +#include + +#include + +namespace picongpu::particles::atomicPhysics::kernel +{ + namespace s_enums = picongpu::particles::atomicPhysics::enums; + + /** update atomic state of ions according to their accepted transition, do nothing if the ion did not accept an + * transition + * + * @note This kernel differs from the RecordChanges kernel by not assuming that all ions accepted a transition. + * @note This kernel is intended for updating the ion atomic state after accepting an instant transition. + */ + template + struct UpdateIonAtomicStateKernel + { + using VectorIdx = pmacc::DataSpace; + + /** call operator + * + * called by UpdateIons atomicPhysics sub-stage for each active processClass + * + * @param worker object containing the device and block information, passed by PMACC_KERNEL call + * @param areaMapping mapping of blockIndex to block superCell index + * @param timeRemainingBox deviceDataBox giving access to the local time remaining of all local super + * cells + * @param ionBox deviceDataBox giving access to the particle frames of all local superCells + * @param atomicStateBox deviceDataBox giving access to atomic state property data + * @param transitionBox deviceDataBox giving access to transition property data, + */ + template< + typename T_Worker, + typename T_AreaMapping, + typename T_LocalTimeRemainingBox, + typename T_IonBox, + typename T_AtomicStateDataBox, + typename T_TransitionDataBox> + HDINLINE void operator()( + T_Worker const& worker, + T_AreaMapping const areaMapping, + T_LocalTimeRemainingBox const timeRemainingBox, + T_IonBox ionBox, + T_AtomicStateDataBox const atomicStateBox, + T_TransitionDataBox const transitionBox) const + { + PMACC_CASSERT(CheckSetOfAtomicDataBoxes:: + transitionDataBoxMatchesProcessClass()); + + VectorIdx const superCellIdx = KernelIndexation::getSuperCellIndex(worker, areaMapping); + VectorIdx const superCellFieldIdx + = KernelIndexation::getSuperCellFieldIndexFromSuperCellIndex(areaMapping, superCellIdx); + + auto const timeRemaining = timeRemainingBox(superCellFieldIdx); + auto forEachLocalIonBoxEntry = pmacc::particles::algorithm::acc::makeForEach(worker, ionBox, superCellIdx); + + bool const superCellAlreadyFinished = (timeRemaining <= 0._X); + bool const superCellHasIons = forEachLocalIonBoxEntry.hasParticles(); + if(superCellAlreadyFinished || !superCellHasIons) + return; + + forEachLocalIonBoxEntry( + [&atomicStateBox, &transitionBox](T_Worker const& worker, auto& ion) + { + bool const accepted = ion[accepted_]; + uint32_t const transitionIndex = ion[transitionIndex_]; + uint8_t const processClass = ion[processClass_]; + + bool const acceptedProcessClassTransition = (processClass != u8(T_ProcessClass)); + + // skip ions without accepted instant transition + if(!accepted || !acceptedProcessClassTransition) + return; + + uint32_t const newAtomicStateCollectionIndex + = transitionBox.upperStateCollectionIndex(transitionIndex); + + picongpu::particles::atomicPhysics::SetAtomicState::op( + atomicStateBox, + ion, + newAtomicStateCollectionIndex); + }); + } + }; +} // namespace picongpu::particles::atomicPhysics::kernel diff --git a/include/picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp b/include/picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp index acb2b5961f..fe80bb43ec 100644 --- a/include/picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp +++ b/include/picongpu/particles/atomicPhysics/rateCalculation/BoundFreeFieldTransitionRates.hpp @@ -29,6 +29,7 @@ #include #include +#include /** @file implements calculation of rates for bound-free field ionization atomic physics transitions * @@ -125,6 +126,8 @@ namespace picongpu::particles::atomicPhysics::rateCalculation * @param Z screenedCharge for ionization electron, e * @param nEff effective principal quantum number, unitless * @param eFieldNorm_AU norm of the E-Field strength, in sim.atomicUnit.eField() + * + * @returns -1. if the result is larger than the maximum representable value of T_ReturnType, unit: 1/unit_time */ template HDINLINE static float_X rateFormula(float_X const Z, float_X const nEff, float_X const eFieldNorm_AU) @@ -138,12 +141,18 @@ namespace picongpu::particles::atomicPhysics::rateCalculation constexpr T_ReturnType pi = pmacc::math::Pi::value; + float_64 subTerm = pmacc::math::cPow(dFromADK, 2u) + * math::exp(float_64(-2._X * ZCubed / (3._X * nEffCubed * eFieldNorm_AU))); + + // explicitly deal with overflows + constexpr auto maxValueT_Return = std::numeric_limits::max(); + if(subTerm > maxValueT_Return) + subTerm = -1.; + // 1/sim.atomicUnit.time() T_ReturnType rateADK_AU = eFieldNorm_AU / (static_cast(8.) * pi * static_cast(Z)) - * static_cast( - pmacc::math::cPow(dFromADK, 2u) - * math::exp(float_64(-2._X * ZCubed / (3._X * nEffCubed * eFieldNorm_AU)))); + * static_cast(subTerm); // factor from averaging over one laser cycle with LINEAR polarization if constexpr( @@ -175,7 +184,7 @@ namespace picongpu::particles::atomicPhysics::rateCalculation * @param atomicStateDataBox access to atomic state property data * @param boundBoundTransitionDataBox access to bound-bound transition data * - * @return unit: 1/picongpu::sim.unit.time() + * @returns -1. if the result is larger than the maximum representable value of T_ReturnType, unit: 1/unit_time */ template< typename T_ReturnType, @@ -221,7 +230,7 @@ namespace picongpu::particles::atomicPhysics::rateCalculation * @param atomicStateDataBox access to atomic state property data * @param boundBoundTransitionDataBox access to bound-bound transition data * - * @return unit: 1/picongpu::sim.unit.time() + * @returns -1. if the result is larger than the maximum representable value of T_ReturnType, unit: 1/unit_time */ template< typename T_ReturnType, diff --git a/include/picongpu/particles/atomicPhysics/stage/CheckForFieldEnergyOverSubscription.hpp b/include/picongpu/particles/atomicPhysics/stage/CheckForFieldEnergyOverSubscription.hpp new file mode 100644 index 0000000000..dc00408b51 --- /dev/null +++ b/include/picongpu/particles/atomicPhysics/stage/CheckForFieldEnergyOverSubscription.hpp @@ -0,0 +1,103 @@ +/* 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 . + */ + +/** @file check for overSubscription of cells and calculate rejectionProbability for each*/ + +#pragma once + +// need picongpu::atomicPhysics::ElectronHistogram from atomicPhysics.param +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldE.hpp" +#include "picongpu/particles/atomicPhysics/electronDistribution/LocalHistogramField.hpp" +#include "picongpu/particles/atomicPhysics/kernel/CheckForFieldEnergyOverSubscription.kernel" +#include "picongpu/particles/atomicPhysics/localHelperFields/FieldEnergyUseCacheField.hpp" +#include "picongpu/particles/atomicPhysics/localHelperFields/RejectionProbabilityCacheField_Cell.hpp" +#include "picongpu/particles/atomicPhysics/localHelperFields/SharedResourcesOverSubscribedField.hpp" +#include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp" + +#include +#include +#include + +#include + +namespace picongpu::particles::atomicPhysics::stage +{ + /** CheckForAndRejectOversubscription atomic physics sub-stage + * + * check each histogram bin for deltaWeight > weight0, if yes mark bin as over subscribed + * + * @tparam T_numberAtomicPhysicsIonSpecies specialization template parameter used to prevent compilation of all + * atomicPhysics kernels if no atomic physics species is present. + */ + template + struct CheckForFieldEnergyOverSubscription + { + //! call of kernel for every superCell + HINLINE void operator()(picongpu::MappingDesc const mappingDesc) const + { + // full local domain, no guards + pmacc::AreaMapping mapper(mappingDesc); + pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); + + auto& timeRemainingField = *dc.get< + picongpu::particles::atomicPhysics::localHelperFields::TimeRemainingField>( + "TimeRemainingField"); + + auto& sharedResourcesOverSubscribedField + = *dc.get>("SharedResourcesOverSubscribedField"); + + auto& rejectionProbabilityCacheField_Cell + = *dc.get>("RejectionProbabilityCacheField_Cell"); + + auto& eField = *dc.get(FieldE::getName()); + + using FieldEnergyUseCacheField + = picongpu::particles::atomicPhysics::localHelperFields::FieldEnergyUseCacheField< + picongpu::MappingDesc>; + auto& fieldEnergyUseCacheField = *dc.get("FieldEnergyUseCacheField"); + + // macro for call of kernel for every superCell, see pull request #4321 + PMACC_LOCKSTEP_KERNEL( + picongpu::particles::atomicPhysics::kernel::CheckForFieldEnergyOverSubscriptionKernel< + T_numberAtomicPhysicsIonSpecies>()) + .template config(mapper.getGridDim())( + mapper, + timeRemainingField.getDeviceDataBox(), + eField.getDeviceDataBox(), + fieldEnergyUseCacheField.getDeviceDataBox(), + sharedResourcesOverSubscribedField.getDeviceDataBox(), + rejectionProbabilityCacheField_Cell.getDeviceDataBox()); + + /// @todo implement photon histogram, Brian Marre, 2023 + } + }; + + //! specialization for no atomicPhysics ion species + template<> + struct CheckForFieldEnergyOverSubscription<0u> + { + //! call of kernel for every superCell + HINLINE void operator()(picongpu::MappingDesc const mappingDesc) const + { + } + }; +} // namespace picongpu::particles::atomicPhysics::stage diff --git a/include/picongpu/particles/atomicPhysics/stage/ChooseInstantTransition.hpp b/include/picongpu/particles/atomicPhysics/stage/ChooseInstantTransition.hpp new file mode 100644 index 0000000000..95ec8a1fa5 --- /dev/null +++ b/include/picongpu/particles/atomicPhysics/stage/ChooseInstantTransition.hpp @@ -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 . + */ + +/** @file ChooseInstantTransition sub-stage of the atomicPhyiscs stage + * + * Selects one transition for all ions with a state loss rate above the instantaneous transition limit + */ + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/fields/FieldE.hpp" +#include "picongpu/particles/atomicPhysics/atomicData/AtomicData.hpp" +#include "picongpu/particles/atomicPhysics/kernel/ChooseInstantNonReversibleTransition.kernel" +#include "picongpu/particles/atomicPhysics/localHelperFields/FoundUnboundIonField.hpp" +#include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp" +#include "picongpu/particles/param.hpp" +#include "picongpu/particles/traits/GetAtomicDataType.hpp" +#include "picongpu/particles/traits/GetNumberAtomicStates.hpp" + +#include +#include +#include + +#include +#include + +namespace picongpu::particles::atomicPhysics::stage +{ + namespace s_enums = picongpu::particles::atomicPhysics::enums; + + /** ChooseInstantTransition atomic physics sub-stage + * + * @tparam T_IonSpecies ion species type + * + * @todo implement reversible instant transitions, Brian Marre, 2024 + * + * @attention assumes that the flag accepted of all macro ions of T_IonSpecies has been reset before calling this + * stage + * @attention assumes that the foundUnboundIonField has been reset before calling this stage + */ + template + struct ChooseInstantTransition + { + // might be alias, from here on out no more + //! resolved type of alias T_IonSpecies + using IonSpecies = pmacc::particles::meta::FindByNameOrType_t; + + // ionization potential depression model to use + using IPDModel = picongpu::atomicPhysics::IPDModel; + + using DistributionFloat = pmacc::random::distributions::Uniform; + using RngFactoryFloat = particles::functor::misc::Rng; + + //! call of kernel for every superCell + HINLINE void operator()( + [[maybe_unused]] picongpu::MappingDesc const mappingDesc, + [[maybe_unused]] uint32_t const currentStep) const + { + using AtomicDataType = typename picongpu::traits::GetAtomicDataType::type; + + if constexpr(AtomicDataType::switchFieldIonization) + { + // full local domain, no guards + pmacc::AreaMapping mapper(mappingDesc); + pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); + + auto& timeRemainingField + = *dc.get>( + "TimeRemainingField"); + auto& foundUnboundIonField = *dc.get< + particles::atomicPhysics::localHelperFields::FoundUnboundIonField>( + "FoundUnboundIonField"); + auto& eField = *dc.get(FieldE::getName()); + auto& ions = *dc.get(IonSpecies::FrameType::getName()); + RngFactoryFloat rngFactoryFloat = RngFactoryFloat{currentStep}; + + constexpr uint32_t numberAtomicStatesOfSpecies + = picongpu::traits::GetNumberAtomicStates::value; + using ChooseInstantNonReversibleTransition = kernel::ChooseInstantNonReversibleTransitionKernel< + IPDModel, + numberAtomicStatesOfSpecies, + AtomicDataType::ADKLaserPolarization>; + + /// @todo implement iteration to capture field strength decrease, Brian Marre, 2024 + auto& atomicData = *dc.get(IonSpecies::FrameType::getName() + "_atomicData"); + IPDModel::template callKernelWithIPDInput< + ChooseInstantNonReversibleTransition, + IonSpecies::FrameType::frameSize>( + dc, + mapper, + rngFactoryFloat, + atomicData.template getChargeStateDataDataBox(), + atomicData.template getAtomicStateDataDataBox(), + atomicData.template getBoundFreeStartIndexBlockDataBox(), + atomicData.template getBoundFreeNumberTransitionsDataBox(), + atomicData + .template getBoundFreeTransitionDataBox(), + timeRemainingField.getDeviceDataBox(), + foundUnboundIonField.getDeviceDataBox(), + eField.getDeviceDataBox(), + ions.getDeviceParticlesBox()); + } + } + }; +} // namespace picongpu::particles::atomicPhysics::stage diff --git a/include/picongpu/particles/atomicPhysics/stage/RecordSuggestedFieldEnergyUse.hpp b/include/picongpu/particles/atomicPhysics/stage/RecordSuggestedFieldEnergyUse.hpp new file mode 100644 index 0000000000..7d331c6511 --- /dev/null +++ b/include/picongpu/particles/atomicPhysics/stage/RecordSuggestedFieldEnergyUse.hpp @@ -0,0 +1,93 @@ +/* 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 . + */ + +//! @file record all accepted transition's suggested changes + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/particles/atomicPhysics/enums/TransitionOrdering.hpp" +#include "picongpu/particles/atomicPhysics/kernel/RecordSuggestedFieldEnergyUse.kernel" +#include "picongpu/particles/atomicPhysics/localHelperFields/FieldEnergyUseCacheField.hpp" +#include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp" +#include "picongpu/particles/param.hpp" + +#include + +#include +#include + +namespace picongpu::particles::atomicPhysics::stage +{ + /** atomicPhysics sub-stage recording for every accepted transition shared physics + * resource usage + * + * for example the histogram in weight usage of a collisional ionization, + * but not the ionization macro electron spawn, since that is not a shared resource. + * + * @attention assumes that the ChooseTransition, ExtractTransitionCollectionIndex + * and AcceptTransitionTest stages have been executed previously in the current + * atomicPhysics time step. + * + * @tparam T_IonSpecies ion species type + */ + template + struct RecordSuggestedFieldEnergyUse + { + // might be alias, from here on out no more + //! resolved type of alias T_IonSpecies + using IonSpecies = pmacc::particles::meta::FindByNameOrType_t; + + //! call of kernel for every superCell + HINLINE void operator()(picongpu::MappingDesc const mappingDesc) const + { + // full local domain, no guards + pmacc::AreaMapping mapper(mappingDesc); + pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); + + using AtomicDataType = typename picongpu::traits::GetAtomicDataType::type; + auto& atomicData = *dc.get(IonSpecies::FrameType::getName() + "_atomicData"); + + auto& timeRemainingField = *dc.get< + picongpu::particles::atomicPhysics::localHelperFields::TimeRemainingField>( + "TimeRemainingField"); + auto& fieldEnergyUseCacheField = *dc.get< + particles::atomicPhysics::localHelperFields::FieldEnergyUseCacheField>( + "FieldEnergyUseCacheField"); + auto& ions = *dc.get(IonSpecies::FrameType::getName()); + + using IPDModel = picongpu::atomicPhysics::IPDModel; + + IPDModel::template callKernelWithIPDInput< + particles::atomicPhysics::kernel::RecordSuggestedFieldEnergyUseKernel, + IonSpecies::FrameType::frameSize>( + dc, + mapper, + atomicData.template getChargeStateDataDataBox(), + atomicData.template getAtomicStateDataDataBox(), + atomicData.template getBoundFreeTransitionDataBox< + false, + picongpu::particles::atomicPhysics::enums::TransitionOrdering::byLowerState>(), + timeRemainingField.getDeviceDataBox(), + fieldEnergyUseCacheField.getDeviceDataBox(), + ions.getDeviceParticlesBox()); + } + }; + +} // namespace picongpu::particles::atomicPhysics::stage diff --git a/include/picongpu/particles/atomicPhysics/stage/UpdateIonAtomicState.hpp b/include/picongpu/particles/atomicPhysics/stage/UpdateIonAtomicState.hpp new file mode 100644 index 0000000000..765d822962 --- /dev/null +++ b/include/picongpu/particles/atomicPhysics/stage/UpdateIonAtomicState.hpp @@ -0,0 +1,84 @@ +/* Copyright 2023-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 . + */ + +//! @file record all ion transitions' delta energy + +#pragma once + +#include "picongpu/defines.hpp" +#include "picongpu/particles/atomicPhysics/electronDistribution/LocalHistogramField.hpp" +#include "picongpu/particles/atomicPhysics/enums/ProcessClass.hpp" +#include "picongpu/particles/atomicPhysics/kernel/UpdateIonAtomicState.kernel" +#include "picongpu/particles/atomicPhysics/localHelperFields/TimeRemainingField.hpp" +#include "picongpu/particles/param.hpp" + +#include + +#include +#include + +namespace picongpu::particles::atomicPhysics::stage +{ + /** atomicPhysics sub-stage updating atomic state according to accepted transitions, does nothing for ions which + * did not accept a transition + * + * @tparam T_IonSpecies ion species type + */ + template + struct UpdateIonAtomicState + { + // might be alias, from here on out no more + //! resolved type of alias T_IonSpecies + using IonSpecies = pmacc::particles::meta::FindByNameOrType_t; + + //! call of kernel for every superCell + HINLINE void operator()(picongpu::MappingDesc const mappingDesc) const + { + using AtomicDataType = typename picongpu::traits::GetAtomicDataType::type; + + if constexpr(AtomicDataType::switchFieldIonization) + { + // full local domain, no guards + pmacc::AreaMapping mapper(mappingDesc); + pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); + + auto& timeRemainingField + = *dc.get>( + "TimeRemainingField"); + auto& ions = *dc.get(IonSpecies::FrameType::getName()); + auto& atomicData = *dc.get(IonSpecies::FrameType::getName() + "_atomicData"); + + namespace s_enums = particles::atomicPhysics::enums; + + using UpdateIonAtomicState_fieldIonization + = picongpu::particles::atomicPhysics::kernel::UpdateIonAtomicStateKernel< + s_enums::ProcessClass::fieldIonization>; + PMACC_LOCKSTEP_KERNEL(UpdateIonAtomicState_fieldIonization()) + .config(mapper.getGridDim(), ions)( + mapper, + timeRemainingField.getDeviceDataBox(), + ions.getDeviceParticlesBox(), + atomicData.template getAtomicStateDataDataBox(), + atomicData.template getBoundFreeTransitionDataBox< + false, + s_enums::TransitionOrdering::byLowerState>()); + } + } + }; +} // namespace picongpu::particles::atomicPhysics::stage diff --git a/include/picongpu/simulation/stage/AtomicPhysics.x.cpp b/include/picongpu/simulation/stage/AtomicPhysics.x.cpp index 77a1d82acf..d8ac5756ef 100644 --- a/include/picongpu/simulation/stage/AtomicPhysics.x.cpp +++ b/include/picongpu/simulation/stage/AtomicPhysics.x.cpp @@ -30,8 +30,10 @@ #include "picongpu/particles/atomicPhysics/param.hpp" #include "picongpu/particles/atomicPhysics/stage/BinElectrons.hpp" #include "picongpu/particles/atomicPhysics/stage/CalculateStepLength.hpp" +#include "picongpu/particles/atomicPhysics/stage/CheckForFieldEnergyOverSubscription.hpp" #include "picongpu/particles/atomicPhysics/stage/CheckForOverSubscription.hpp" #include "picongpu/particles/atomicPhysics/stage/CheckPresence.hpp" +#include "picongpu/particles/atomicPhysics/stage/ChooseInstantTransition.hpp" #include "picongpu/particles/atomicPhysics/stage/ChooseTransition.hpp" #include "picongpu/particles/atomicPhysics/stage/ChooseTransitionGroup.hpp" #include "picongpu/particles/atomicPhysics/stage/DecelerateElectrons.hpp" @@ -40,6 +42,7 @@ #include "picongpu/particles/atomicPhysics/stage/LoadAtomicInputData.hpp" #include "picongpu/particles/atomicPhysics/stage/RecordChanges.hpp" #include "picongpu/particles/atomicPhysics/stage/RecordSuggestedChanges.hpp" +#include "picongpu/particles/atomicPhysics/stage/RecordSuggestedFieldEnergyUse.hpp" #include "picongpu/particles/atomicPhysics/stage/ResetAcceptedStatus.hpp" #include "picongpu/particles/atomicPhysics/stage/ResetRateCache.hpp" #include "picongpu/particles/atomicPhysics/stage/ResetSharedResources.hpp" @@ -47,6 +50,7 @@ #include "picongpu/particles/atomicPhysics/stage/RollForOverSubscription.hpp" #include "picongpu/particles/atomicPhysics/stage/SpawnIonizationElectrons.hpp" #include "picongpu/particles/atomicPhysics/stage/UpdateElectricField.hpp" +#include "picongpu/particles/atomicPhysics/stage/UpdateIonAtomicState.hpp" #include "picongpu/particles/atomicPhysics/stage/UpdateTimeRemaining.hpp" #include "picongpu/particles/filter/filter.hpp" @@ -77,7 +81,8 @@ namespace picongpu::simulation::stage { SubStep, ChooseTransition, - RejectOverSubscription + RejectOverSubscription, + ApplyInstantTransitions }; } // namespace enums @@ -108,7 +113,6 @@ namespace picongpu::simulation::stage uint32_t T_numberAtomicPhysicsIonSpecies> struct AtomicPhysics { - private: // linearized dataBox of SuperCellField template using LinearizedBox = DataBoxDim1Access; @@ -132,6 +136,7 @@ namespace picongpu::simulation::stage using IPDIonSpecies = MakeSeq_t; //!@} + private: //! debug print to console //!@{ @@ -230,14 +235,6 @@ namespace picongpu::simulation::stage } //!@} - //! set timeRemaining to PIC-time step - HINLINE static void setTimeRemaining() - { - pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); - auto& localTimeRemainingField = *dc.get("TimeRemainingField"); - localTimeRemainingField.getDeviceBuffer().setValue(picongpu::sim.pic.getDt()); // sim.unit.time() - } - //! reset the histogram on device side HINLINE static void resetElectronEnergyHistogram() { @@ -249,14 +246,20 @@ namespace picongpu::simulation::stage } //! reset foundUnboundIonField on device side - HINLINE static void resetFoundUnboundIon() + template + HINLINE static void resetFoundUnboundIon(T_FoundUnboundIonField& foundUnboundIonField) { - pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); - auto& foundUnboundIonField = *dc.get("FoundUnboundIonField"); foundUnboundIonField.getDeviceBuffer().setValue(0._X); }; - HINLINE static void resetAcceptStatus(picongpu::MappingDesc const& mappingDesc) + //! reset SharedResourcesOverSubscribedField on device side + HINLINE static void resetSharedResourceOverSubscribed(picongpu::MappingDesc const& mappingDesc) + { + picongpu::particles::atomicPhysics::stage::ResetSharedResources{}( + mappingDesc); + }; + + HINLINE static void resetAcceptedStatus(picongpu::MappingDesc const& mappingDesc) { // particle[accepted_] = false, in each macro ion using ForEachIonSpeciesResetAcceptedStatus = pmacc::meta::ForEach< @@ -265,6 +268,30 @@ namespace picongpu::simulation::stage ForEachIonSpeciesResetAcceptedStatus{}(mappingDesc); } + //! reset each superCell's time step + HINLINE static void resetTimeStep(picongpu::MappingDesc const& mappingDesc) + { + // timeStep = timeRemaining + picongpu::particles::atomicPhysics::stage::ResetTimeStepField()( + mappingDesc); + } + + //! reset each superCell's rate cache + HINLINE static void resetRateCache() + { + using ForEachIonSpeciesResetRateCache = pmacc::meta:: + ForEach>; + ForEachIonSpeciesResetRateCache{}(); + } + + //! set timeRemaining to PIC-time step + HINLINE static void setTimeRemaining() + { + pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); + auto& localTimeRemainingField = *dc.get("TimeRemainingField"); + localTimeRemainingField.getDeviceBuffer().setValue(picongpu::sim.pic.getDt()); // sim.unit.time() + } + HINLINE static void debugForceConstantElectronTemperature([[maybe_unused]] uint32_t const currentStep) { if constexpr(debug::scFlyComparison::FORCE_CONSTANT_ELECTRON_TEMPERATURE) @@ -296,22 +323,6 @@ namespace picongpu::simulation::stage currentStep); } - //! reset each superCell's time step - HINLINE static void resetTimeStep(picongpu::MappingDesc const& mappingDesc) - { - // timeStep = timeRemaining - picongpu::particles::atomicPhysics::stage::ResetTimeStepField()( - mappingDesc); - } - - //! reset each superCell's rate cache - HINLINE static void resetRateCache() - { - using ForEachIonSpeciesResetRateCache = pmacc::meta:: - ForEach>; - ForEachIonSpeciesResetRateCache{}(); - } - //! check which atomic states are actually present in each superCell HINLINE static void checkPresence(picongpu::MappingDesc const& mappingDesc) { @@ -344,9 +355,9 @@ namespace picongpu::simulation::stage ForEachIonSpeciesCalculateStepLength{}(mappingDesc); } + //! randomly roll transition for each not yet accepted macro ion HINLINE static void chooseTransition(picongpu::MappingDesc const& mappingDesc, uint32_t const currentStep) { - // randomly roll transition for each not yet accepted macro ion using ForEachIonSpeciesChooseTransitionGroup = pmacc::meta::ForEach< AtomicPhysicsIonSpecies, particles::atomicPhysics::stage::ChooseTransitionGroup>; @@ -358,18 +369,36 @@ namespace picongpu::simulation::stage ForEachIonSpeciesChooseTransition{}(mappingDesc, currentStep); } + //! randomly roll transition for each macro ion with an instant transition + HINLINE static void chooseInstantTransition( + picongpu::MappingDesc const& mappingDesc, + uint32_t const currentStep) + { + using ForEachIonSpeciesChooseInstantTransition = pmacc::meta::ForEach< + AtomicPhysicsIonSpecies, + picongpu::particles::atomicPhysics::stage::ChooseInstantTransition>; + ForEachIonSpeciesChooseInstantTransition{}(mappingDesc, currentStep); + } + //! record all shared resources usage by accepted transitions HINLINE static void recordSuggestedChanges(picongpu::MappingDesc const& mappingDesc) { - picongpu::particles::atomicPhysics::stage::ResetSharedResources{}( - mappingDesc); using ForEachIonSpeciesRecordSuggestedChanges = pmacc::meta::ForEach< AtomicPhysicsIonSpecies, particles::atomicPhysics::stage::RecordSuggestedChanges>; ForEachIonSpeciesRecordSuggestedChanges{}(mappingDesc); } - // check if an electron histogram bin () is over subscription --> superCellOversubScriptionField + //! record all shared resources usage by accepted transitions + HINLINE static void recordSuggestedFieldEnergyUse(picongpu::MappingDesc const& mappingDesc) + { + using ForEachIonSpeciesRecordSuggestedFieldEnergyUse = pmacc::meta::ForEach< + AtomicPhysicsIonSpecies, + particles::atomicPhysics::stage::RecordSuggestedFieldEnergyUse>; + ForEachIonSpeciesRecordSuggestedFieldEnergyUse{}(mappingDesc); + } + + // check if an electron histogram bin or cell's electric field is over subscribed template< enums::Loop T_Loop, typename T_SuperCellSharedResourcesOverSubscriptionField, @@ -379,32 +408,74 @@ namespace picongpu::simulation::stage T_SuperCellSharedResourcesOverSubscriptionField& perSuperCellSharedResourcesOverSubscriptionField, T_DeviceReduce& deviceReduce) { - DataSpace const fieldGridLayoutOverSubscription - = perSuperCellSharedResourcesOverSubscriptionField.getGridLayout().sizeWithoutGuardND(); - + resetSharedResourceOverSubscribed(mappingDesc); picongpu::particles::atomicPhysics::stage::CheckForOverSubscription{}( mappingDesc); + DataSpace const fieldGridLayoutOverSubscription + = perSuperCellSharedResourcesOverSubscriptionField.getGridLayout().sizeWithoutGuardND(); auto linearizedOverSubscribedBox = LinearizedBox( perSuperCellSharedResourcesOverSubscriptionField.getDeviceDataBox(), fieldGridLayoutOverSubscription); - bool isOverSubscribed = static_cast(deviceReduce( + bool const isOverSubscribed = static_cast(deviceReduce( pmacc::math::operation::Or(), linearizedOverSubscribedBox, fieldGridLayoutOverSubscription.productOfComponents())); // debug only - std::string message; if constexpr(debugPrintActive()) { + std::string message; message += "[histogram oversubscribed?]: "; message += (isOverSubscribed ? "true" : "false"); + + printHistogramToConsole(mappingDesc, message); + printFieldEnergyUseCacheToConsole(mappingDesc); + printOverSubscriptionFieldToConsole(mappingDesc); + printRejectionProbabilityCacheToConsole(mappingDesc); + } + + // check whether a least one histogram is oversubscribed + return isOverSubscribed; + } + + // check if a cell's electric field is over subscribed + template< + enums::Loop T_Loop, + typename T_SuperCellSharedResourcesOverSubscriptionField, + typename T_DeviceReduce> + HINLINE static bool isElectricFieldOverSubscribed( + picongpu::MappingDesc const& mappingDesc, + T_SuperCellSharedResourcesOverSubscriptionField& perSuperCellSharedResourcesOverSubscriptionField, + T_DeviceReduce& deviceReduce) + { + resetSharedResourceOverSubscribed(mappingDesc); + picongpu::particles::atomicPhysics::stage::CheckForFieldEnergyOverSubscription< + T_numberAtomicPhysicsIonSpecies>{}(mappingDesc); + + DataSpace const fieldGridLayoutOverSubscription + = perSuperCellSharedResourcesOverSubscriptionField.getGridLayout().sizeWithoutGuardND(); + auto linearizedOverSubscribedBox = LinearizedBox( + perSuperCellSharedResourcesOverSubscriptionField.getDeviceDataBox(), + fieldGridLayoutOverSubscription); + + bool const isOverSubscribed = static_cast(deviceReduce( + pmacc::math::operation::Or(), + linearizedOverSubscribedBox, + fieldGridLayoutOverSubscription.productOfComponents())); + + // debug only + if constexpr(debugPrintActive()) + { + std::string message = "[histogram oversubscribed?]: "; + message += (isOverSubscribed ? "true" : "false"); + + printHistogramToConsole(mappingDesc, message); + printFieldEnergyUseCacheToConsole(mappingDesc); + printOverSubscriptionFieldToConsole(mappingDesc); + printRejectionProbabilityCacheToConsole(mappingDesc); } - printHistogramToConsole(mappingDesc, message); - printFieldEnergyUseCacheToConsole(mappingDesc); - printOverSubscriptionFieldToConsole(mappingDesc); - printRejectionProbabilityCacheToConsole(mappingDesc); // check whether a least one histogram is oversubscribed return isOverSubscribed; @@ -423,7 +494,8 @@ namespace picongpu::simulation::stage /** update atomic state and accumulate delta energy for delta energy histogram * * @note may already update the atomic state since the following kernels DecelerateElectrons and - * SpawnIonizationElectrons only use the transitionIndex particle attribute */ + * SpawnIonizationElectrons only use the transitionIndex particle attribute + */ HINLINE static void recordChanges(picongpu::MappingDesc const& mappingDesc) { using ForEachIonSpeciesRecordChanges = pmacc::meta:: @@ -431,6 +503,25 @@ namespace picongpu::simulation::stage ForEachIonSpeciesRecordChanges{}(mappingDesc); } + //! update atomic state of all ions having accepted an instant transition + HINLINE static void updateIonAtomicState(picongpu::MappingDesc const& mappingDesc) + { + using ForEachIonSpeciesUpdateAtomicState = pmacc::meta::ForEach< + AtomicPhysicsIonSpecies, + particles::atomicPhysics::stage::UpdateIonAtomicState>; + ForEachIonSpeciesUpdateAtomicState{}(mappingDesc); + } + + HINLINE static void spawnIonizationElectrons( + picongpu::MappingDesc const& mappingDesc, + uint32_t const currentStep) + { + using ForEachIonSpeciesSpawnIonizationElectrons = pmacc::meta::ForEach< + AtomicPhysicsIonSpecies, + particles::atomicPhysics::stage::SpawnIonizationElectrons>; + ForEachIonSpeciesSpawnIonizationElectrons{}(mappingDesc, currentStep); + } + HINLINE static void updateElectrons(picongpu::MappingDesc const& mappingDesc, uint32_t const currentStep) { /** @note DecelerateElectrons must be called before SpawnIonizationElectrons such that we only @@ -440,10 +531,7 @@ namespace picongpu::simulation::stage particles::atomicPhysics::stage::DecelerateElectrons>; ForEachElectronSpeciesDecelerateElectrons{}(mappingDesc); - using ForEachIonSpeciesSpawnIonizationElectrons = pmacc::meta::ForEach< - AtomicPhysicsIonSpecies, - particles::atomicPhysics::stage::SpawnIonizationElectrons>; - ForEachIonSpeciesSpawnIonizationElectrons{}(mappingDesc, currentStep); + spawnIonizationElectrons(mappingDesc, currentStep); } HINLINE static void updateElectricField(picongpu::MappingDesc const& mappingDesc) @@ -453,7 +541,7 @@ namespace picongpu::simulation::stage } template - HINLINE static void doIPDIonization( + HINLINE static void applyIPDIonization( picongpu::MappingDesc const& mappingDesc, uint32_t const currentStep, T_DeviceReduce& deviceReduce) @@ -468,7 +556,7 @@ namespace picongpu::simulation::stage bool foundUnbound = true; do { - resetFoundUnboundIon(); + resetFoundUnboundIon(foundUnboundIonField); calculateIPDInput(mappingDesc, currentStep); picongpu::atomicPhysics::IPDModel::template applyIPDIonization( mappingDesc, @@ -486,6 +574,63 @@ namespace picongpu::simulation::stage while(foundUnbound); } + template + HINLINE static void applyInstantTransitions( + picongpu::MappingDesc const& mappingDesc, + uint32_t const currentStep, + T_DeviceReduce& deviceLocalReduce) + { + pmacc::DataConnector& dc = pmacc::Environment<>::get().DataConnector(); + + /// @todo pass instead of duplicating code from applyIPDIonization?, Brian Marre, 2024 + auto& foundUnboundIonField = *dc.get("FoundUnboundIonField"); + DataSpace const fieldGridLayoutFoundUnbound + = foundUnboundIonField.getGridLayout().sizeWithoutGuardND(); + + auto& perSuperCellSharedResourcesOverSubscribedField + = *dc.get("SharedResourcesOverSubscribedField"); + + // instant Transition loop, ends when no ion in state with instant transition anymore + bool foundInstantTransitionIon; + do + { + resetFoundUnboundIon(foundUnboundIonField); + chooseInstantTransition(mappingDesc, currentStep); + recordSuggestedFieldEnergyUse(mappingDesc); + + bool isFieldOverSubscribed = isElectricFieldOverSubscribed( + mappingDesc, + perSuperCellSharedResourcesOverSubscribedField, + deviceLocalReduce); + + while(isFieldOverSubscribed) + { + // at least one cell's field energy over-subscribed + randomlyRejectTransitionFromOverSubscribedResources(mappingDesc, currentStep); + recordSuggestedFieldEnergyUse(mappingDesc); + + isFieldOverSubscribed = isElectricFieldOverSubscribed( + mappingDesc, + perSuperCellSharedResourcesOverSubscribedField, + deviceLocalReduce); + } // end remove over subscription loop + + updateIonAtomicState(mappingDesc); + spawnIonizationElectrons(mappingDesc, currentStep); + + updateElectricField(mappingDesc); + + auto linearizedFoundUnboundIonBox = LinearizedBox( + foundUnboundIonField.getDeviceDataBox(), + fieldGridLayoutFoundUnbound); + foundInstantTransitionIon = static_cast(deviceLocalReduce( + pmacc::math::operation::Or(), + linearizedFoundUnboundIonBox, + fieldGridLayoutFoundUnbound.productOfComponents())); + } // end instant transition loop + while(foundInstantTransitionIon); + } + HINLINE static void updateTimeRemaining(picongpu::MappingDesc const& mappingDesc) { // timeRemaining -= timeStep @@ -539,10 +684,12 @@ namespace picongpu::simulation::stage bool isSubSteppingComplete = false; while(!isSubSteppingComplete) { - resetAcceptStatus(mappingDesc); + resetAcceptedStatus(mappingDesc); resetElectronEnergyHistogram(); debugForceConstantElectronTemperature(currentStep); - doIPDIonization(mappingDesc, currentStep, deviceLocalReduce); + applyIPDIonization(mappingDesc, currentStep, deviceLocalReduce); + applyInstantTransitions(mappingDesc, currentStep, deviceLocalReduce); + resetAcceptedStatus(mappingDesc); binElectronsToEnergyHistogram(mappingDesc); resetTimeStep(mappingDesc); resetRateCache(); @@ -591,7 +738,7 @@ namespace picongpu::simulation::stage } // end atomicPhysics sub-stepping loop // ensure no unbound states are visible to the rest of the loop - doIPDIonization(mappingDesc, currentStep, deviceLocalReduce); + applyIPDIonization(mappingDesc, currentStep, deviceLocalReduce); } }; diff --git a/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics.param b/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics.param index a94701a36e..42614d3c03 100644 --- a/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics.param +++ b/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics.param @@ -45,16 +45,25 @@ namespace picongpu::atomicPhysics MaxEnergyParam>; //} - // atomicPhysics rate solver settings - + //! 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.1_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 diff --git a/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics_Debug.param b/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics_Debug.param index 5e76b6f641..0dc2832316 100644 --- a/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics_Debug.param +++ b/share/picongpu/examples/AtomicPhysics/include/picongpu/param/atomicPhysics_Debug.param @@ -139,8 +139,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