Skip to content

Commit

Permalink
Merge pull request #5218 from BrianMarre/topic-instantFieldIonization
Browse files Browse the repository at this point in the history
instant field ionization transitions
  • Loading branch information
psychocoderHPC authored Dec 9, 2024
2 parents db6f900 + 5b182ae commit 63999f7
Show file tree
Hide file tree
Showing 22 changed files with 1,458 additions and 235 deletions.
31 changes: 20 additions & 11 deletions include/picongpu/param/atomicPhysics.param
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,6 @@

#pragma once

#include "picongpu/particles/atomicPhysics/ExponentialApproximationProbability.hpp"
#include "picongpu/particles/atomicPhysics/LinearApproximationProbability.hpp"
#include "picongpu/particles/atomicPhysics/atomicData/AtomicData.def"
#include "picongpu/particles/atomicPhysics/electronDistribution/LogSpaceHistogram.hpp"
#include "picongpu/particles/atomicPhysics/enums/ADKLaserPolarization.hpp"
Expand All @@ -34,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
Expand Down
8 changes: 8 additions & 0 deletions include/picongpu/param/atomicPhysics_Debug.param
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

This file was deleted.

Original file line number Diff line number Diff line change
@@ -0,0 +1,45 @@
/* Copyright 2024 Brian Marre
*
* This file is part of PIConGPU.
*
* PIConGPU is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* PIConGPU is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PIConGPU.
* If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once

// need unit system and normalization
#include "picongpu/defines.hpp"
#include "picongpu/particles/atomicPhysics/param.hpp"

namespace picongpu::particles::atomicPhysics
{
struct InstantTransitionRateLimit
{
/** get maximum of total state loss rate for inclusion in the time dependent rate equation solver
*
* @tparam T_ReturnType type and precision to use in the result
*/
template<typename T_ReturnType>
static constexpr T_ReturnType get()
{
using picongpu::atomicPhysics::RateSolverParam;

// unit: unitless * unitless / unit_time = 1/unit_time
return static_cast<T_ReturnType>(
RateSolverParam::timeStepAlpha * float_X(RateSolverParam::maximumNumberSubStepsPerPICTimeStep))
/ picongpu::sim.pic.getDt<T_ReturnType>();
}
};
} // namespace picongpu::particles::atomicPhysics

This file was deleted.

Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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{});
});
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
/* Copyright 2024 Brian Marre
*
* This file is part of PIConGPU.
*
* PIConGPU is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* PIConGPU is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PIConGPU.
* If not, see <http://www.gnu.org/licenses/>.
*/

#pragma once

#include "picongpu/defines.hpp"
#include "picongpu/particles/atomicPhysics/ConvertEnum.hpp"
#include "picongpu/particles/atomicPhysics/KernelIndexation.hpp"
#include "picongpu/particles/atomicPhysics/kernel/CalculateRejectionProbability.hpp"
#include "picongpu/particles/atomicPhysics/localHelperFields/RejectionProbabilityCache_Cell.hpp"

#include <pmacc/lockstep/ForEach.hpp>

#include <cstdint>

namespace picongpu::particles::atomicPhysics::kernel
{
/** check all cells for overSubscription
*
* checks each superCell cell for fieldEnergyUse > fieldEnergy,
* if yes mark cell as oversubscribed and stores rejection probability in rejection probability cache
*
* @tparam T_numberAtomicPhysicsIonSpecies specialization template parameter used to prevent compilation of all
* atomicPhysics kernels if no atomic physics species is present.
*/
template<uint32_t T_numberAtomicPhysicsIonSpecies>
struct CheckForFieldEnergyOverSubscriptionKernel
{
using VectorIdx = DataSpace<picongpu::simDim>;
/** call operator
*
* called by CheckForOverSubscription atomic physics sub-stage
*
* @param worker object containing the device and block information, passed by PMACC_KERNEL call
* @param areMapping mapping of blockIndex to block superCell index
* @param rngFactory factory for uniformly distributed random number generator
* @param eFieldBox deviceDataBox giving access to eField values for all local superCells
* @param fieldEnergyUseCacheBox deviceDataBox giving access to the field energy use cache for each local
* superCell
* @param sharedResourcesOverSubscribedBox deviceDataBox giving access to the local shared resources over
* subscription switch for each local superCell
* @param rejectionProbabilityCacheCellBox deviceDataBox giving access to localRejectionProbabilityCache for
* all cells of all local superCells
*/
template<
typename T_Worker,
typename T_AreaMapping,
typename T_TimeRemainingBox,
typename T_EFieldBox,
typename T_FieldEnergyUseCacheBox,
typename T_SharedRessourcesOverSubscribedBox,
typename T_RejectionProbabilityCacheCellDataBox>
HDINLINE void operator()(
T_Worker const& worker,
T_AreaMapping const areaMapping,
T_TimeRemainingBox const timeRemainingBox,
T_EFieldBox const eFieldBox,
T_FieldEnergyUseCacheBox const fieldEnergyUseCacheBox,
T_SharedRessourcesOverSubscribedBox sharedResourcesOverSubscribedBox,
T_RejectionProbabilityCacheCellDataBox rejectionProbabilityCacheCellBox) const
{
auto const superCellIdx = KernelIndexation::getSuperCellIndex(worker, areaMapping);
auto const superCellFieldIdx
= KernelIndexation::getSuperCellFieldIndexFromSuperCellIndex(areaMapping, superCellIdx);

auto const timeRemaining = timeRemainingBox(superCellFieldIdx);
if(timeRemaining <= 0._X)
return;

using FieldEnergyCache = typename T_FieldEnergyUseCacheBox::ValueType;
using RejectionProbabilityCache
= localHelperFields::RejectionProbabilityCache_Cell<FieldEnergyCache::numberCells>;
RejectionProbabilityCache& rejectionProbabilityCache = rejectionProbabilityCacheCellBox(superCellFieldIdx);

FieldEnergyCache const& eFieldEnergyUseCache = fieldEnergyUseCacheBox(superCellFieldIdx);

VectorIdx const superCellCellOffset = superCellIdx * picongpu::SuperCellSize::toRT();
bool sharedResourcesOverSubscribed = false;

auto forEachCell = pmacc::lockstep::makeForEach<FieldEnergyCache::numberCells, T_Worker>(worker);
forEachCell(
[&worker,
&superCellCellOffset,
&eFieldBox,
&eFieldEnergyUseCache,
&rejectionProbabilityCache,
&sharedResourcesOverSubscribed](uint32_t const linearCellIndex)
{
if(CalculateRejectionProbability::ofCell(
linearCellIndex,
superCellCellOffset,
eFieldBox,
eFieldEnergyUseCache,
rejectionProbabilityCache))
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
Loading

0 comments on commit 63999f7

Please sign in to comment.