diff --git a/src/picongpu/include/particles/Particles.hpp b/src/picongpu/include/particles/Particles.hpp index cbd16c3991..d2d3357686 100644 --- a/src/picongpu/include/particles/Particles.hpp +++ b/src/picongpu/include/particles/Particles.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2013-2014 Axel Huebl, Heiko Burau, Rene Widera, Felix Schmitt, + * Copyright 2013-2015 Axel Huebl, Heiko Burau, Rene Widera, Felix Schmitt, * Marco Garten * * This file is part of PIConGPU. @@ -72,10 +72,6 @@ class Particles : public ParticlesBase, publ void syncToDevice(); - //Ionization - template - void ionize(T_Elec electrons, uint32_t currentStep); - private: SimulationDataId datasetID; GridLayout gridLayout; diff --git a/src/picongpu/include/particles/Particles.tpp b/src/picongpu/include/particles/Particles.tpp index 8652c1e05a..0397d5744f 100644 --- a/src/picongpu/include/particles/Particles.tpp +++ b/src/picongpu/include/particles/Particles.tpp @@ -1,5 +1,5 @@ /** - * Copyright 2013-2014 Axel Huebl, Heiko Burau, Rene Widera, Richard Pausch, Felix Schmitt + * Copyright 2013-2015 Axel Huebl, Heiko Burau, Rene Widera, Richard Pausch, Felix Schmitt * * This file is part of PIConGPU. * @@ -30,9 +30,6 @@ #include "particles/Particles.kernel" -// Ionization -#include "particles/ionization/ionization.hpp" - #include "dataManagement/DataConnector.hpp" #include "mappings/kernel/AreaMapping.hpp" diff --git a/src/picongpu/include/particles/ParticlesFunctors.hpp b/src/picongpu/include/particles/ParticlesFunctors.hpp index e2fdf2137b..52ba0b9ea3 100644 --- a/src/picongpu/include/particles/ParticlesFunctors.hpp +++ b/src/picongpu/include/particles/ParticlesFunctors.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Rene Widera, Marco Garten + * Copyright 2014-2015 Rene Widera, Marco Garten * * This file is part of PIConGPU. * @@ -155,29 +155,74 @@ struct CallUpdate } }; -/* Tests if species can be ionized and calls the function to do that */ +/** \struct CallIonization + * + * \brief Tests if species can be ionized and calls the kernel to do that + * + * \tparam T_SpeciesName name of particle species that is checked for ionization + */ template struct CallIonization { typedef T_SpeciesName SpeciesName; typedef typename SpeciesName::type SpeciesType; - + typedef typename SpeciesType::FrameType FrameType; + /* SelectIonizer will be either the specified one or fallback: None */ typedef typename GetIonizer::type SelectIonizer; - /* describes the instance of CallIonization */ - template + /** Functor implementation + * + * \tparam T_StorageStuple contains info about the particle species + * \tparam T_CellDescription contains the number of blocks and blocksize + * that is later passed to the kernel + * \param tuple An n-tuple containing the type-info of multiple particle species + * \param cellDesc points to logical block information like dimension and cell sizes + * \param currentStep The current time step + */ + template HINLINE void operator()( T_StorageTuple& tuple, + T_CellDescription* cellDesc, const uint32_t currentStep ) const { - - /* alias for pointer on source species */ - PMACC_AUTO(speciesPtr, tuple[SpeciesName()]); - /* instance of particle ionizer that was flagged in speciesDefinition.param */ - SelectIonizer myIonizer; - myIonizer(*speciesPtr, tuple, currentStep); - + /* only if an ionizer has been specified, this is executed */ + typedef typename HasFlag >::type hasIonizer; + if (hasIonizer::value) + { + + /* define the type of the species to be created + * from inside the ionization model specialization + */ + typedef typename SelectIonizer::DestSpecies DestSpecies; + /* alias for pointer on source species */ + PMACC_AUTO(srcSpeciesPtr, tuple[SpeciesName()]); + /* alias for pointer on destination species */ + PMACC_AUTO(electronsPtr, tuple[typename MakeIdentifier::type()]); + + /* 3-dim vector : number of threads to be started in every dimension */ + dim3 block( MappingDesc::SuperCellSize::toRT().toDim3() ); + + /** kernelIonizeParticles + * \brief calls the ionization model and handles that electrons are created correctly + * while cycling through the particle frames + * + * kernel call : instead of name<<>> (args, ...) + * "blocks" will be calculated from "this->cellDescription" and "CORE + BORDER" + * "threads" is calculated from the previously defined vector "block" + */ + __picKernelArea( particles::ionization::kernelIonizeParticles, *cellDesc, CORE + BORDER ) + (block) + ( srcSpeciesPtr->getDeviceParticlesBox( ), + electronsPtr->getDeviceParticlesBox( ), + SelectIonizer(currentStep) + ); + /* fill the gaps in the created species' particle frames to ensure that only + * the last frame is not completely filled but every other before is full + */ + electronsPtr->fillAllGaps(); + + } } }; // struct CallIonization diff --git a/src/picongpu/include/particles/ionization/None/AlgorithmNone.hpp b/src/picongpu/include/particles/ionization/None/AlgorithmNone.hpp index 9559eff15f..f9c34d6f96 100644 --- a/src/picongpu/include/particles/ionization/None/AlgorithmNone.hpp +++ b/src/picongpu/include/particles/ionization/None/AlgorithmNone.hpp @@ -22,9 +22,13 @@ #include "types.h" -/** IONIZATION ALGORITHM +/** \file AlgorithNone.hpp + * + * IONIZATION ALGORITHM for the model None + * * - implements the calculation of ionization probability and changes charge states - * - is called with the IONIZATION MODEL, specifically by setting the flag in @see speciesDefinition.param */ + * - is called with the IONIZATION MODEL, specifically by setting the flag in @see speciesDefinition.param + */ namespace picongpu { @@ -34,11 +38,14 @@ namespace ionization { /** \struct AlgorithmNone - * \brief ionization algorithm that does nothing */ + * + * \brief ionization algorithm that does nothing + */ struct AlgorithmNone { /** Functor implementation + * * \tparam EType type of electric field * \tparam BType type of magnetic field * \tparam ParticleType type of particle to be ionized diff --git a/src/picongpu/include/particles/ionization/None/None.def b/src/picongpu/include/particles/ionization/None/None.def index 77ddfae3a7..05245d6b16 100644 --- a/src/picongpu/include/particles/ionization/None/None.def +++ b/src/picongpu/include/particles/ionization/None/None.def @@ -20,6 +20,8 @@ #pragma once +#include "types.h" + namespace picongpu { namespace particles @@ -27,9 +29,35 @@ namespace particles namespace ionization { + /** \struct None_Impl + * + * \brief Empty ionization model [that is used for all species that are not ionized] + * + * \tparam T_DestSpecies electron species to be created + * \tparam T_SrcSpecies particle species that is ionized + * default is boost::mpl placeholder because specialization + * cannot be known in list of particle species' flags + * \see speciesDefinition.param + */ + template + struct None_Impl; + /** \struct None - * \brief Fallback for all species that cannot/should not be ionized */ - struct None; + * + * \brief Fallback for all species that cannot/should not be ionized + * + * \tparam T_DestSpecies electron species to be created + * + * wrapper class, + * needed because the SrcSpecies cannot be known during the + * first specialization of the ionization model in the particle definition + * \see speciesDefinition.param + */ + template + struct None + { + typedef None_Impl type; + }; } // namespace ionization } // namespace particles diff --git a/src/picongpu/include/particles/ionization/None/None.hpp b/src/picongpu/include/particles/ionization/None/None_Impl.hpp similarity index 61% rename from src/picongpu/include/particles/ionization/None/None.hpp rename to src/picongpu/include/particles/ionization/None/None_Impl.hpp index a03c4a8e62..5a4f6f7777 100644 --- a/src/picongpu/include/particles/ionization/None/None.hpp +++ b/src/picongpu/include/particles/ionization/None/None_Impl.hpp @@ -20,6 +20,7 @@ #pragma once +#include "types.h" #include "particles/ionization/None/None.def" #include "particles/ionization/None/AlgorithmNone.hpp" @@ -29,17 +30,35 @@ namespace particles { namespace ionization { - + /* fallback for all species that cannot/should not be ionized */ - struct None + template + struct None_Impl { - template - void operator()(T_SrcSpecies& src, T_ParticleStorage& pst, const uint32_t currentStep) - { - // Do nothing - } + typedef T_DestSpecies DestSpecies; + typedef T_SrcSpecies SrcSpecies; + + private: + + public: + + None_Impl(const uint32_t currentStep) + { + + } + + DINLINE void init(const DataSpace& blockCell, const int& linearThreadIdx, const DataSpace& totalCellOffset) const + { + + } + + template + DINLINE void operator()(FrameType& ionFrame, int localIdx, unsigned int& newMacroElectrons) const + { + + } }; - + } // namespace ionization } // namespace particles } // namespace picongpu diff --git a/src/picongpu/include/particles/ionization/byField/BSI/AlgorithmBSI.hpp b/src/picongpu/include/particles/ionization/byField/BSI/AlgorithmBSI.hpp index 68b3e5e46a..bbf4903eee 100644 --- a/src/picongpu/include/particles/ionization/byField/BSI/AlgorithmBSI.hpp +++ b/src/picongpu/include/particles/ionization/byField/BSI/AlgorithmBSI.hpp @@ -26,10 +26,14 @@ #include "traits/attribute/GetChargeState.hpp" #include "algorithms/math/floatMath/floatingPoint.tpp" -/** IONIZATION ALGORITHM +/** \file AlgorithmBSI.hpp + * + * IONIZATION ALGORITHM for the BSI model + * * - implements the calculation of ionization probability and changes charge states * by decreasing the number of bound electrons - * - is called with the IONIZATION MODEL, specifically by setting the flag in @see speciesDefinition.param */ + * - is called with the IONIZATION MODEL, specifically by setting the flag in @see speciesDefinition.param + */ namespace picongpu { @@ -39,11 +43,14 @@ namespace ionization { /** \struct AlgorithmBSI - * \brief calculation for the Barrier Suppression Ionization model */ + * + * \brief calculation for the Barrier Suppression Ionization model + */ struct AlgorithmBSI { - + /** Functor implementation + * * \tparam EType type of electric field * \tparam BType type of magnetic field * \tparam ParticleType type of particle to be ionized @@ -59,9 +66,9 @@ namespace ionization const float_X protonNumber = GetAtomicNumbers::type::numberOfProtons; float_X chargeState = attribute::getChargeState(parentIon); - + uint32_t cs = math::float2int_rd(chargeState); - + /* ionization condition */ if (math::abs(eField)*UNIT_EFIELD >= SI::IONIZATION_EFIELD[cs] && chargeState < protonNumber) { diff --git a/src/picongpu/include/particles/ionization/byField/BSI/BSI.def b/src/picongpu/include/particles/ionization/byField/BSI/BSI.def index 134658ad27..2924785770 100644 --- a/src/picongpu/include/particles/ionization/byField/BSI/BSI.def +++ b/src/picongpu/include/particles/ionization/byField/BSI/BSI.def @@ -20,14 +20,29 @@ #pragma once +#include "types.h" + namespace picongpu { namespace particles { namespace ionization { + /** \struct BSI_Impl + * + * \brief Barrier Suppression Ionization - Implementation + * + * \tparam T_DestSpecies electron species to be created + * \tparam T_SrcSpecies particle species that is ionized + * default is boost::mpl placeholder because specialization + * cannot be known in list of particle species' flags + * \see speciesDefinition.param + */ + template + struct BSI_Impl; /** \struct BSI + * * \brief Barrier Suppression Ionization * * - takes the ionization energies of the various charge states of ions @@ -35,9 +50,18 @@ namespace ionization * - if the field strength is locally exceeded: increase the charge state * - see for example: Delone, N. B.; Krainov, V. P. (1998). "Tunneling and barrier-suppression ionization of atoms and ions in a laser radiation field" * - * \tparam T_DestSpecies electron species to be created */ + * \tparam T_DestSpecies electron species to be created + * + * wrapper class, + * needed because the SrcSpecies cannot be known during the + * first specialization of the ionization model in the particle definition + * \see speciesDefinition.param + */ template - struct BSI; + struct BSI + { + typedef BSI_Impl type; + }; } // namespace ionization } // namespace particles diff --git a/src/picongpu/include/particles/ionization/byField/BSI/BSI.hpp b/src/picongpu/include/particles/ionization/byField/BSI/BSI.hpp deleted file mode 100644 index c16f65832b..0000000000 --- a/src/picongpu/include/particles/ionization/byField/BSI/BSI.hpp +++ /dev/null @@ -1,64 +0,0 @@ -/** - * Copyright 2015 Marco Garten - * - * 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 "particles/ionization/byField/BSI/BSI.def" -#include "particles/ionization/byField/BSI/AlgorithmBSI.hpp" - -#include "compileTime/conversion/TypeToPointerPair.hpp" - -namespace picongpu -{ -namespace particles -{ -namespace ionization -{ - - /** \struct BSI - * \brief Barrier Suppression Ionization - */ - template - struct BSI - { - - /* define ionization ALGORITHM for ionization MODEL */ - typedef particles::ionization::AlgorithmBSI IonizationAlgorithm; - - /** Functor implementation - * \tparam T_SrcSpecies particle species to be ionized - * \tparam T_ParticleStorage \see picongpu/src/libPMacc/include/math/MapTuple.hpp - * - * \param currentStep current time step */ - template - void operator()(T_SrcSpecies& src, T_ParticleStorage& pst, const uint32_t currentStep) - { - - /* alias for electron species to be created */ - PMACC_AUTO(electrons,pst[typename MakeIdentifier::type()]); - /* call the ionization function */ - src.ionize(electrons, currentStep); - } - - }; - -} // namespace ionization -} // namespace particles -} // namespace picongpu diff --git a/src/picongpu/include/particles/ionization/byField/BSI/BSI_Impl.hpp b/src/picongpu/include/particles/ionization/byField/BSI/BSI_Impl.hpp new file mode 100644 index 0000000000..3d8f3df7d0 --- /dev/null +++ b/src/picongpu/include/particles/ionization/byField/BSI/BSI_Impl.hpp @@ -0,0 +1,193 @@ +/** + * Copyright 2015 Marco Garten + * + * 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 "types.h" +#include "math/vector/Size_t.hpp" +#include "simulation_defines.hpp" +#include "traits/Resolve.hpp" +#include "mappings/kernel/AreaMapping.hpp" + +#include "fields/FieldB.hpp" +#include "fields/FieldE.hpp" + +#include "particles/ionization/byField/BSI/BSI.def" +#include "particles/ionization/byField/BSI/AlgorithmBSI.hpp" +#include "particles/ionization/ionization.hpp" + +#include "compileTime/conversion/TypeToPointerPair.hpp" +#include "memory/boxes/DataBox.hpp" + +#include "particles/ParticlesFunctors.hpp" + +namespace picongpu +{ +namespace particles +{ +namespace ionization +{ + + /** \struct BSI_Impl + * \brief Barrier Suppression Ionization - Implementation + * + * \tparam T_DestSpecies electron species to be created + * \tparam T_SrcSpecies particle species that is ionized + */ + template + struct BSI_Impl + { + + typedef T_DestSpecies DestSpecies; + typedef T_SrcSpecies SrcSpecies; + + typedef typename SrcSpecies::FrameType FrameType; + + /* specify field to particle interpolation scheme */ + typedef typename PMacc::traits::Resolve< + typename GetFlagType >::type + >::type Field2ParticleInterpolation; + + /* margins around the supercell for the interpolation of the field on the cells */ + typedef typename GetMargin::LowerMargin LowerMargin; + typedef typename GetMargin::UpperMargin UpperMargin; + + /* relevant area of a block */ + typedef SuperCellDescription< + typename MappingDesc::SuperCellSize, + LowerMargin, + UpperMargin + > BlockArea; + + BlockArea BlockDescription; + + private: + + /* define ionization ALGORITHM (calculation) for ionization MODEL */ + typedef particles::ionization::AlgorithmBSI IonizationAlgorithm; + + typedef MappingDesc::SuperCellSize TVec; + + typedef FieldE::ValueType ValueType_E; + typedef FieldB::ValueType ValueType_B; + /* global memory EM-field device databoxes */ + FieldE::DataBoxType eBox; + FieldB::DataBoxType bBox; + /* shared memory EM-field device databoxes */ + PMACC_ALIGN(cachedE, DataBox >); + PMACC_ALIGN(cachedB, DataBox >); + + public: + /* host constructor */ + BSI_Impl(const uint32_t currentStep) + { + DataConnector &dc = Environment<>::get().DataConnector(); + /* initialize pointers on host-side E-(B-)field databoxes */ + FieldE* fieldE = &(dc.getData (FieldE::getName(), true)); + FieldB* fieldB = &(dc.getData (FieldB::getName(), true)); + /* initialize device-side E-(B-)field databoxes */ + eBox = fieldE->getDeviceDataBox(); + bBox = fieldB->getDeviceDataBox(); + + } + + /** Initialization function on device + * + * \brief Cache EM-fields on device + * and initialize possible prerequisites for ionization, like e.g. random number generator. + * + * This function will be called inline on the device which must happen BEFORE threads diverge + * during loop execution. The reason for this is the `__syncthreads()` call which is necessary after + * initializing the E-/B-field shared boxes in shared memory. + */ + DINLINE void init(const DataSpace& blockCell, const int& linearThreadIdx, const DataSpace& totalCellOffset) + { + + /* caching of E and B fields */ + cachedB = CachedBox::create < 0, ValueType_B > (BlockArea()); + cachedE = CachedBox::create < 1, ValueType_E > (BlockArea()); + /* wait for shared memory to be initialized */ + __syncthreads(); + /* instance of nvidia assignment operator */ + nvidia::functors::Assign assign; + /* copy fields from global to shared */ + PMACC_AUTO(fieldBBlock, bBox.shift(blockCell)); + ThreadCollective collective(linearThreadIdx); + collective( + assign, + cachedB, + fieldBBlock + ); + /* copy fields from global to shared */ + PMACC_AUTO(fieldEBlock, eBox.shift(blockCell)); + collective( + assign, + cachedE, + fieldEBlock + ); + } + + /** Functor implementation + * + * \param ionFrame reference to frame of the to-be-ionized particles + * \param localIdx local (linear) index in super cell / frame + * \param newMacroElectrons reference to variable for each thread that stores the number + * of macro electrons to be created during the current time step + */ + DINLINE void operator()(FrameType& ionFrame, int localIdx, unsigned int& newMacroElectrons) + { + + /* type of PIC-scheme cell */ + typedef typename fieldSolver::NumericalCellType NumericalCellType; + + /* alias for the single macro-particle */ + PMACC_AUTO(particle,ionFrame[localIdx]); + /* particle position, used for field-to-particle interpolation */ + floatD_X pos = particle[position_]; + const int particleCellIdx = particle[localCellIdx_]; + /* multi-dim coordinate of the local cell inside the super cell */ + DataSpace localCell(DataSpaceOperations::template map (particleCellIdx)); + /* interpolation of E- */ + ValueType_E eField = Field2ParticleInterpolation() + (cachedE.shift(localCell).toCursor(), pos, NumericalCellType::getEFieldPosition()); + /* and B-field on the particle position */ + ValueType_B bField = Field2ParticleInterpolation() + (cachedB.shift(localCell).toCursor(), pos, NumericalCellType::getBFieldPosition()); + + /* define number of bound macro electrons before ionization */ + float_X prevBoundElectrons = particle[boundElectrons_]; + + /* this is the point where actual ionization takes place */ + IonizationAlgorithm ionizeAlgo; + ionizeAlgo( + bField, eField, + particle + ); + + /* determine number of new macro electrons to be created */ + newMacroElectrons = prevBoundElectrons - particle[boundElectrons_]; + + } + + }; + +} // namespace ionization +} // namespace particles +} // namespace picongpu diff --git a/src/picongpu/include/particles/ionization/byField/ionizers.def b/src/picongpu/include/particles/ionization/byField/ionizers.def index a00de6b84b..c699639696 100644 --- a/src/picongpu/include/particles/ionization/byField/ionizers.def +++ b/src/picongpu/include/particles/ionization/byField/ionizers.def @@ -20,7 +20,10 @@ #pragma once -/* These includes contain forward declarations of < Ionization Models > */ +/* \file ionizers.def + * + * These includes contain forward declarations of < Ionization Models > + */ #include "particles/ionization/byField/BSI/BSI.def" #include "particles/ionization/None/None.def" diff --git a/src/picongpu/include/particles/ionization/byField/ionizers.hpp b/src/picongpu/include/particles/ionization/byField/ionizers.hpp index a94e613535..cb75405f30 100644 --- a/src/picongpu/include/particles/ionization/byField/ionizers.hpp +++ b/src/picongpu/include/particles/ionization/byField/ionizers.hpp @@ -20,9 +20,12 @@ #pragma once -/* Includes containing definition of < Ionization Models > +/** \file ionizers.hpp + * + * Includes containing definition of < Ionization Models > * which itself each include their own < Ionization Algorithm > - * that implements what the model actually DOES */ + * that implements what the model actually DOES + */ -#include "particles/ionization/byField/BSI/BSI.hpp" -#include "particles/ionization/None/None.hpp" +#include "particles/ionization/byField/BSI/BSI_Impl.hpp" +#include "particles/ionization/None/None_Impl.hpp" diff --git a/src/picongpu/include/particles/ionization/ionization.hpp b/src/picongpu/include/particles/ionization/ionization.hpp index fc6419c5a1..1af458f7ca 100644 --- a/src/picongpu/include/particles/ionization/ionization.hpp +++ b/src/picongpu/include/particles/ionization/ionization.hpp @@ -1,6 +1,6 @@ /** - * Copyright 2014 Marco Garten, Axel Huebl, Heiko Burau, Rene Widera, - * Richard Pausch, Felix Schmitt + * Copyright 2014-2015 Marco Garten, Axel Huebl, Heiko Burau, Rene Widera, + * Richard Pausch, Felix Schmitt * * This file is part of PIConGPU. * @@ -38,110 +38,58 @@ namespace picongpu { - -using namespace PMacc; - -/** \struct IonizeParticlesPerFrame - * \brief gathers fields for particle interpolation and gives them to the ionization model - * - * \tparam IonizeAlgo ionization algorithm - * \tparam TVec dimensions of the target: - * e.g. size of the super cell - * \tparam T_Field2ParticleInterpolation type of field to particle interpolation - * \tparam NumericalCellType type of cell with respect to field discretization: - e.g. Yee cell - */ -template -struct IonizeParticlesPerFrame +namespace particles +{ +namespace ionization { - /** Functor implementation - * - * \tparam FrameType frame type of particles that get ionized - * \tparam T_BBox type of B-field box - * \tparam T_EBox type of E-field box - * - * \param ionFrame (here address of: ) frame of the to-be-ionized particles - * \param localIdx local (linear) index in super cell / frame - * \param bBox B-field box instance - * \param ebox E-field box instance - * \param newMacroElectrons (here address of: ) variable for each thread that stores the number - * of macro electrons to be created during the current time step - */ - template - DINLINE void operator()(FrameType& ionFrame, int localIdx, T_BBox& bBox, T_EBox& eBox, unsigned int& newMacroElectrons) - { - - typedef TVec Block; - /** type for field to particle interpolation */ - typedef T_Field2ParticleInterpolation Field2ParticleInterpolation; - - typedef typename T_BBox::ValueType BType; - typedef typename T_EBox::ValueType EType; - - PMACC_AUTO(particle,ionFrame[localIdx]); - - floatD_X pos = particle[position_]; - const int particleCellIdx = particle[localCellIdx_]; - - DataSpace localCell(DataSpaceOperations::template map (particleCellIdx)); - - - EType eField = Field2ParticleInterpolation() - (eBox.shift(localCell).toCursor(), pos, NumericalCellType::getEFieldPosition()); - - BType bField =Field2ParticleInterpolation() - (bBox.shift(localCell).toCursor(), pos, NumericalCellType::getBFieldPosition()); - - /* define number of bound macro electrons before ionization */ - float_X prevBoundElectrons = particle[boundElectrons_]; - - /* this is the point where actual ionization takes place */ - IonizeAlgo ionizeAlgo; - ionizeAlgo( - bField, eField, - particle - ); - - /* determine number of new macro electrons to be created */ - newMacroElectrons = prevBoundElectrons - particle[boundElectrons_]; - /*calculate one dimensional cell index*/ - - } -}; // struct IonizeParticlesPerFrame +using namespace PMacc; /** kernelIonizeParticles * \brief main kernel for ionization * * - maps the frame dimensions and gathers the particle boxes - * - caches the E- and B-fields * - contains / calls the ionization algorithm * - calls the electron creation functors * - * \tparam BlockDescription_ container for information about block / super cell dimensions * \tparam ParBoxIons container of the ions * \tparam ParBoxElectrons container of the electrons - * \tparam BBox b-field box class - * \tparam EBox e-field box class * \tparam Mapping class containing methods for acquiring info from the block - * \tparam FrameIonizer \see IonizeParticlesPerFrame + * \tparam FrameIonizer \see e.g. BSI_Impl in BSI_Impl.hpp + * instance of the ionization model functor */ -template +template __global__ void kernelIonizeParticles(ParBoxIons ionBox, ParBoxElectrons electronBox, - EBox fieldE, - BBox fieldB, FrameIonizer frameIonizer, Mapping mapper) { - + /* "particle box" : container/iterator where the particles live in - * and where one can get the frame in a super cell from */ + * and where one can get the frame in a super cell from + */ typedef typename ParBoxElectrons::FrameType ELECTRONFRAME; typedef typename ParBoxIons::FrameType IONFRAME; - /* for not mixing assign up with the nVidia functor assign */ + + /* specify field to particle interpolation scheme */ + typedef typename PMacc::traits::Resolve< + typename GetFlagType >::type + >::type InterpolationScheme; + + /* margins around the supercell for the interpolation of the field on the cells */ + typedef typename GetMargin::LowerMargin LowerMargin; + typedef typename GetMargin::UpperMargin UpperMargin; + + /* relevant area of a block */ + typedef SuperCellDescription< + typename MappingDesc::SuperCellSize, + LowerMargin, + UpperMargin + > BlockDescription_; + + /* for not mixing operations::assign up with the nvidia functor assign */ namespace partOp = PMacc::particles::operations; - + /* definitions for domain variables, like indices of blocks and threads */ typedef typename BlockDescription_::SuperCellSize SuperCellSize; /* "offset" 3D distance to origin in units of super cells */ @@ -154,7 +102,11 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, /* "offset" from origin of the grid in unit of cells */ const DataSpace blockCell = block * SuperCellSize::toRT(); - + + /* subtract guarding cells to only have the simulation volume */ + const DataSpace localCellIndex = (block * SuperCellSize::toRT() + threadIndex) - mapper.getGuardingSuperCells() * SuperCellSize::toRT(); + + /* typedef for the functor that writes new macro electrons into electron frames during runtime */ typedef typename particles::ionization::WriteElectronIntoFrame WriteElectronIntoFrame; __syncthreads(); @@ -167,7 +119,8 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, __syncthreads(); /*wait that all shared memory is initialized*/ /* find last frame in super cell - * define maxParticlesInFrame as the maximum frame size */ + * define maxParticlesInFrame as the maximum frame size + */ if (linearThreadIdx == 0) { ionFrame = &(ionBox.getLastFrame(block, isValid)); @@ -178,46 +131,29 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, if (!isValid) return; //end kernel if we have no frames - /* caching of E and B fields */ - PMACC_AUTO(cachedB, CachedBox::create < 0, typename BBox::ValueType > (BlockDescription_())); - PMACC_AUTO(cachedE, CachedBox::create < 1, typename EBox::ValueType > (BlockDescription_())); + /* caching of E- and B- fields and initialization of random generator if needed */ + frameIonizer.init(blockCell, linearThreadIdx, localCellIndex); - __syncthreads(); - - nvidia::functors::Assign assign; - - PMACC_AUTO(fieldBBlock, fieldB.shift(blockCell)); - ThreadCollective collective(linearThreadIdx); - collective( - assign, - cachedB, - fieldBBlock - ); - - PMACC_AUTO(fieldEBlock, fieldE.shift(blockCell)); - collective( - assign, - cachedE, - fieldEBlock - ); - /* Declare counter in shared memory that will later tell the current fill level or - * occupation of the newly created target electron frames. */ + * occupation of the newly created target electron frames. + */ __shared__ int newFrameFillLvl; - - __syncthreads(); /*wait that all shared memory is initialized*/ - + + __syncthreads(); /*wait until all shared memory is initialized*/ + /* Declare local variable oldFrameFillLvl for each thread */ int oldFrameFillLvl; - + /* Initialize local (register) counter for each thread - * - describes how many new macro electrons should be created */ + * - describes how many new macro electrons should be created + */ unsigned int newMacroElectrons = 0; - + /* Declare local electron ID - * - describes at which position in the new frame the new electron is to be created */ + * - describes at which position in the new frame the new electron is to be created + */ int electronId; - + /* Master initializes the frame fill level with 0 */ if (linearThreadIdx == 0) { @@ -225,12 +161,13 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, electronFrame = NULL; } __syncthreads(); - + /* move over source species frames and call frameIonizer * frames are worked on in backwards order to avoid asking if there is another frame * --> performance * Because all frames are completely filled except the last and apart from that last frame - * one wants to make sure that all threads are working and every frame is worked on. */ + * one wants to make sure that all threads are working and every frame is worked on. + */ while (isValid) { /* casting uint8_t multiMask to boolean */ @@ -239,11 +176,12 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, /* < IONIZATION and change of charge states > * if the threads contain particles, the frameIonizer can ionize them - * if they are non-particles their inner ionization counter remains at 0 */ + * if they are non-particles their inner ionization counter remains at 0 + */ if (isParticle) /* ionization based on ionization model - this actually increases charge states*/ - frameIonizer(*ionFrame, linearThreadIdx, cachedB, cachedE, newMacroElectrons); - + frameIonizer(*ionFrame, linearThreadIdx, newMacroElectrons); + __syncthreads(); /* always true while-loop over all particles inside source frame until each thread breaks out individually * @@ -251,14 +189,16 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, * The question might arise what happens if more electrons are created than would fit into two frames. * Well, multi-ionization during a time step is accounted for. The number of new electrons is * determined inside the outer loop over the valid frames while in the inner loop each thread can create only ONE - * new macro electron. But the loop repeats until each thread has created all the electrons needed in the time step. */ + * new macro electron. But the loop repeats until each thread has created all the electrons needed in the time step. + */ while (true) { /* < INIT > * - electronId is initialized as -1 (meaning: invalid) * - (local) oldFrameFillLvl set equal to (shared) newFrameFillLvl for each thread * --> each thread remembers the old "counter" - * - then sync */ + * - then sync + */ electronId = -1; oldFrameFillLvl = newFrameFillLvl; __syncthreads(); @@ -266,21 +206,23 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, * - if a thread wants to create electrons in each cycle it can do that only once * and before that it atomically adds to the shared counter and uses the current * value as electronId in the new frame - * - then sync */ + * - then sync + */ if (newMacroElectrons > 0) electronId = atomicAdd(&newFrameFillLvl, 1); - + __syncthreads(); /* < EXIT? > * - if the counter hasn't changed all threads break out of the loop */ if (oldFrameFillLvl == newFrameFillLvl) break; - + __syncthreads(); /* < FIRST NEW FRAME > * - if there is no frame, yet, the master will create a new target electron frame * and attach it to the back of the frame list - * - sync all threads again for them to know which frame to use */ + * - sync all threads again for them to know which frame to use + */ if (linearThreadIdx == 0) { if (electronFrame == NULL) @@ -293,19 +235,21 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, /* < CREATE 1 > * - all electrons fitting into the current frame are created there * - internal ionization counter is decremented by 1 - * - sync */ + * - sync + */ if ((0 <= electronId) && (electronId < maxParticlesInFrame)) { /* each thread makes the attributes of its ion accessible */ PMACC_AUTO(parentIon,((*ionFrame)[linearThreadIdx])); /* each thread initializes an electron if one should be created */ PMACC_AUTO(targetElectronFull,((*electronFrame)[electronId])); - + /* create an electron in the new electron frame: - * - see particles/ionization/ionizationMethods.hpp */ + * - see particles/ionization/ionizationMethods.hpp + */ WriteElectronIntoFrame writeElectron; writeElectron(parentIon,targetElectronFull); - + newMacroElectrons -= 1; } __syncthreads(); @@ -313,7 +257,8 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, * - if the shared counter is larger than the frame size a new electron frame is reserved * and attached to the back of the frame list * - then the shared counter is set back by one frame size - * - sync so that every thread knows about the new frame */ + * - sync so that every thread knows about the new frame + */ if (linearThreadIdx == 0) { if (newFrameFillLvl >= maxParticlesInFrame) @@ -328,7 +273,8 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, * - if the EID is larger than the frame size * - the EID is set back by one frame size * - the thread writes an electron to the new frame - * - the internal counter is decremented by 1 */ + * - the internal counter is decremented by 1 + */ if (electronId >= maxParticlesInFrame) { electronId -= maxParticlesInFrame; @@ -337,12 +283,13 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, PMACC_AUTO(parentIon,((*ionFrame)[linearThreadIdx])); /* each thread initializes an electron if one should be produced */ PMACC_AUTO(targetElectronFull,((*electronFrame)[electronId])); - + /* create an electron in the new electron frame: - * - see particles/ionization/ionizationMethods.hpp */ + * - see particles/ionization/ionizationMethods.hpp + */ WriteElectronIntoFrame writeElectron; writeElectron(parentIon,targetElectronFull); - + newMacroElectrons -= 1; } __syncthreads(); @@ -358,60 +305,6 @@ __global__ void kernelIonizeParticles(ParBoxIons ionBox, } } // void kernelIonizeParticles -/** ionize - * \brief ionization function currently being a member function of the source species - * - * \tparam T_ParticleDescription container of particle source species information - * \tparam T_Elec type of species to be created during ionization - */ -template< typename T_ParticleDescription > -template< typename T_Elec > -void Particles::ionize( T_Elec electrons, uint32_t ) -{ - - /* get the alias for the used ionizer (ionization model) and specify the ionization algorithm */ - typedef typename GetFlagType >::type IonizerAlias; - typedef typename PMacc::traits::Resolve::type::IonizationAlgorithm IonizationAlgorithm; - - typedef typename PMacc::traits::Resolve< - typename GetFlagType >::type - >::type InterpolationScheme; - - /* margins around the supercell for the interpolation of the field on the cells */ - typedef typename GetMargin::LowerMargin LowerMargin; - typedef typename GetMargin::UpperMargin UpperMargin; - - /* frame ionizer that moves over the frames with the particles and executes an operation on them */ - typedef IonizeParticlesPerFrame FrameIonizer; - - /* relevant area of a block */ - typedef SuperCellDescription< - typename MappingDesc::SuperCellSize, - LowerMargin, - UpperMargin - > BlockArea; - - /* 3-dim vector : number of threads to be started in every dimension */ - dim3 block( MappingDesc::SuperCellSize::toRT().toDim3() ); - - /* kernel call : instead of name<<>> (args, ...) - "blocks" will be calculated from "this->cellDescription" and "CORE + BORDER" - "threads" is calculated from the previously defined vector "block" */ - __picKernelArea( kernelIonizeParticles, this->cellDescription, CORE + BORDER ) - (block) - ( this->getDeviceParticlesBox( ), - electrons->getDeviceParticlesBox(), - this->fieldE->getDeviceDataBox( ), - this->fieldB->getDeviceDataBox( ), - FrameIonizer( ) - ); - - /* fill the gaps in the created species' particle frames to ensure that only - * the last frame is not completely filled but every other before is full */ - electrons->fillAllGaps(); - -} // Particles::ionize - +} // namespace ionization +} // namespace particles } // namespace picongpu diff --git a/src/picongpu/include/particles/ionization/ionizationMethods.hpp b/src/picongpu/include/particles/ionization/ionizationMethods.hpp index 70bab0cf14..3c01287fa1 100644 --- a/src/picongpu/include/particles/ionization/ionizationMethods.hpp +++ b/src/picongpu/include/particles/ionization/ionizationMethods.hpp @@ -35,10 +35,10 @@ namespace ionization { using namespace PMacc; - + /* operations on particles */ namespace partOp = PMacc::particles::operations; - + /** \struct WriteElectronIntoFrame * * \brief functor that fills an electron frame entry with details about the created particle @@ -63,7 +63,8 @@ namespace ionization * - multiMask: reading from global memory takes longer than just setting it again explicitly * - momentum: because the electron would get a higher energy because of the ion mass * - boundElectrons: because species other than ions or atoms do not have them - * (gets AUTOMATICALLY deselected because electrons do not have this attribute) */ + * (gets AUTOMATICALLY deselected because electrons do not have this attribute) + */ PMACC_AUTO(targetElectronClone, partOp::deselect >(childElectron)); partOp::assign(targetElectronClone, parentIon); @@ -74,16 +75,16 @@ namespace ionization float3_X electronMomentum (parentIon[momentum_]*(massElectron/massIon)); childElectron[momentum_] = electronMomentum; - + /* conservation of momentum * \todo add conservation of mass */ parentIon[momentum_] -= electronMomentum; } }; - + } // namespace ionization } // namespace particles - + } // namespace picongpu diff --git a/src/picongpu/include/particles/traits/GetIonizer.hpp b/src/picongpu/include/particles/traits/GetIonizer.hpp index 45cb32c474..e2e9309503 100644 --- a/src/picongpu/include/particles/traits/GetIonizer.hpp +++ b/src/picongpu/include/particles/traits/GetIonizer.hpp @@ -1,5 +1,5 @@ /** - * Copyright 2014 Marco Garten + * Copyright 2014-2015 Marco Garten * * This file is part of PIConGPU. * @@ -20,6 +20,7 @@ #pragma once +#include "types.h" #include "simulation_defines.hpp" #include "particles/memory/frames/Frame.hpp" #include "traits/GetFlagType.hpp" @@ -27,25 +28,37 @@ #include "simulation_defines/param/speciesDefinition.param" #include "simulation_defines/unitless/speciesDefinition.unitless" +#include "particles/ionization/byField/ionizers.def" +#include "particles/ionization/byField/ionizers.hpp" + namespace picongpu { -template +template struct GetIonizer { - - typedef typename T_Species::FrameType FrameType; + + typedef T_SpeciesType SpeciesType; + typedef typename SpeciesType::FrameType FrameType; typedef typename HasFlag >::type hasIonizer; - + /* The following line only fetches the alias */ typedef typename GetFlagType >::type FoundIonizerAlias; + /* This now resolves the alias into the actual object type */ typedef typename PMacc::traits::Resolve::type FoundIonizer; - /* if no ionizer was defined we use IonizerNone as fallback */ - typedef typename bmpl::if_::type type; - -}; -}// namespace picongpu + /* This specifies the source species as the second template parameter of the ionization model */ + typedef typename bmpl::if_< + hasIonizer, + FoundIonizer, + particles::ionization::None + >::type UserIonizer; + + /* specializes the designated ionization model with the particle species it is called upon */ + typedef typename bmpl::apply1::type type; +}; // struct GetIonizer + +}// namespace picongpu diff --git a/src/picongpu/include/simulationControl/MySimulation.hpp b/src/picongpu/include/simulationControl/MySimulation.hpp index 676d4315e5..6afeef591c 100644 --- a/src/picongpu/include/simulationControl/MySimulation.hpp +++ b/src/picongpu/include/simulationControl/MySimulation.hpp @@ -384,9 +384,11 @@ class MySimulation : public SimulationHelper { namespace nvfct = PMacc::nvidia::functors; - /* Ionization */ + /* Initialize ionization routine for each species + * - valid species will be ionized + * - invalid species (e.g. electrons): fallback */ ForEach, MakeIdentifier > particleIonization; - particleIonization(forward(particleStorage), currentStep); + particleIonization(forward(particleStorage), cellDescription, currentStep); EventTask initEvent = __getTransactionEvent(); EventTask updateEvent;