From 30f6fc38fd2962cd00c85356602427cdce981f73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Filip=20Opto=C5=82owicz?= Date: Thu, 6 Jun 2024 13:21:52 +0200 Subject: [PATCH 1/2] Example change - doesn't work yet. Just need to create another pull request --- .../synchrotron/AlgorithmSynchrotron.hpp | 668 ++++++++++-------- 1 file changed, 366 insertions(+), 302 deletions(-) diff --git a/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp b/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp index a25905f52b..054297fe9a 100644 --- a/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp +++ b/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp @@ -71,358 +71,422 @@ namespace picongpu F = math::pow(F, 1._X - f) * math::pow(Ftemp, f); } - /** Functor implementation - * @tparam EType type of electric field - * @tparam BType type of magnetic field - * @tparam ParticleType type of electron particle - * - * @param bField magnetic field value at t=0 - * @param eField electric field value at t=0 - * @param parentElectron instance of electron with position at t=0 and momentum at t=-1/2 - * @param randNr random number, equally distributed in range [0.:1.0] - * - * @return photon energy with respect to electron energy if photon is created 0 otherwise - */ - + //! Retrieves the interpolation parameters used for calculating zq and its index + HDINLINE constexpr auto getInterpolationParams() + { + // return a tuple with the following values: + // - minZqExponent: minimum zq exponent + // - maxZqExponent: maximum zq exponent + // - numberTableEntries: number of entries in the interpolation table + // - stepWidthLogarithmicScale: step width for logarithmic scale + return std::make_tuple( + params::InterpolationParams::minZqExponent, + params::InterpolationParams::maxZqExponent, + static_cast(params::InterpolationParams::numberTableEntries - 1), + (static_cast( + params::InterpolationParams::maxZqExponent - params::InterpolationParams::minZqExponent) + / static_cast(params::InterpolationParams::numberTableEntries - 1))); + } + //! Calculates the value of delta from the random number randNr1 + HDINLINE float_X calculateDelta(float_64 randNr1) + { + float_X r1r1 = randNr1 * randNr1; + return r1r1 * randNr1; + } + //! Calculates the effective field HeffValue based on particle's velocity, bField, and eField template - HDINLINE float_X operator()( - const T_Worker& worker, - const BType bField, - const EType eField, - ParticleType& parentElectron, - float_64 randNr1, - float_64 randNr2, - GridBuffer::DataBoxType F1F2DeviceBuff, - GridBuffer::DataBoxType failedRequirementQBuff) const + HDINLINE float_X + calculateHeff(const ParticleType& parentElectron, const BType& bField, const EType& eField) { - constexpr float_64 minZqExponent = params::InterpolationParams::minZqExponent; - constexpr float_64 maxZqExponent = params::InterpolationParams::maxZqExponent; - constexpr float_X interpolationPoints - = static_cast(params::InterpolationParams::numberTableEntries - 1); - constexpr float_X stepWidthLogatihmicScale - = (static_cast(maxZqExponent - minZqExponent) / interpolationPoints); - - /** zq = 2/3 * chi^(-1) * delta / (1 - delta) ; - * chi - ratio of the typical photon energy to the electron kinetic energy: - * if larger than 1 the radiation is in the quantum regime - * delta - the ratio of the photon energy to the electron energy - * see Extended particle-in-cell schemes for physics in ultrastrong laser fields: Review and - * developments, A. Gonoskov et al. */ - - //! delta = randNr1^3 but we use randNr1^2 later so we calculate it here - float_X const r1r1 = randNr1 * randNr1; - float_X const delta = r1r1 * randNr1; - - //! calculate Heff - float_X const mass = attribute::getMass(parentElectron[weighting_], parentElectron); - float3_X const vel = Velocity()(parentElectron[momentum_], mass); - - float_X const Vmag = pmacc::math::l2norm(vel); - float3_X const crossVB = pmacc::math::cross(vel, bField); - //! Velocity normalized - float3_X const Vnorm = vel / Vmag; - float_X const dotVnormE = pmacc::math::dot(Vnorm, eField); - float3_X const eFieldPlusCrossVB = eField + crossVB; - + float_X mass = attribute::getMass(parentElectron[weighting_], parentElectron); + float3_X vel = Velocity()(parentElectron[momentum_], mass); + float_X Vmag = pmacc::math::l2norm(vel); + float3_X crossVB = pmacc::math::cross(vel, bField); + float3_X Vnorm = vel / Vmag; + float_X dotVnormE = pmacc::math::dot(Vnorm, eField); + float3_X eFieldPlusCrossVB = eField + crossVB; float_X HeffValue = pmacc::math::dot(eFieldPlusCrossVB, eFieldPlusCrossVB) - dotVnormE * dotVnormE; - /* this check is important to avoid numerical rounding errors that result in NaN - * when taking a sqrt - */ + if(HeffValue <= 0) { return 0; } - HeffValue = math::sqrt(HeffValue); - - - //! Calculate chi - float_X const gamma = Gamma()(parentElectron[momentum_], mass); - float_X const inverse_E_Schwinger = (-ELECTRON_CHARGE * HBAR) + return math::sqrt(HeffValue); + } + //! Calculates the quantum parameter chi + template + HDINLINE float_X calculateChi(const ParticleType& parentElectron, float_X HeffValue) + { + float_X mass = attribute::getMass(parentElectron[weighting_], parentElectron); + float_X gamma = Gamma()(parentElectron[momentum_], mass); + float_X inverse_E_Schwinger = (-ELECTRON_CHARGE * HBAR) / (ELECTRON_MASS * ELECTRON_MASS * SPEED_OF_LIGHT * SPEED_OF_LIGHT * SPEED_OF_LIGHT); - float_X const chi = HeffValue * gamma * inverse_E_Schwinger; - - //! zq - const float_64 oneMinusDeltaOverDelta = (1. - delta) / delta; - const float_X zq = 2._X / (3._X * chi) / oneMinusDeltaOverDelta; - - //! zq convert to index and F1 and F2 - const float_X zqExponent = math::log2(zq); - const int16_t index - = static_cast((zqExponent - minZqExponent) / stepWidthLogatihmicScale); - + return HeffValue * gamma * inverse_E_Schwinger; + } + //! Calculates the zq value and its corresponding index for interpolation + HDINLINE std::tuple calculateZq( + float_X chi, + float_X delta, + float_64 minZqExponent, + float_X stepWidthLogarithmicScale) + { + float_64 oneMinusDeltaOverDelta = (1. - delta) / delta; + float_X zq = 2._X / (3._X * chi) / oneMinusDeltaOverDelta; + float_X zqExponent = math::log2(zq); + int16_t index = static_cast((zqExponent - minZqExponent) / stepWidthLogarithmicScale); + return std::make_tuple(zq, index); + } + //! Interpolates the values of F1 and F2 based on zq and the precomputed interpolation parameters + HDINLINE std::tuple interpolateF1F2( + int16_t index, + float_X zq, + float_X minZqExponent, + float_X interpolationPoints, + float_X stepWidthLogarithmicScale, + GridBuffer::DataBoxType F1F2DeviceBuff) + { float_X F1 = 0._X; float_X F2 = 0._X; float_X Ftemp = 0._X; if(index >= 0 && index < interpolationPoints - 2) { - float_X zq1 = math::pow(2, minZqExponent + stepWidthLogatihmicScale * (index)); - float_X zq2 = math::pow(2, minZqExponent + stepWidthLogatihmicScale * (index + 1)); + float_X zq1 = math::pow(2, minZqExponent + stepWidthLogarithmicScale * (index)); + float_X zq2 = math::pow(2, minZqExponent + stepWidthLogarithmicScale * (index + 1)); float_X f = (zq - zq1) / (zq2 - zq1); - /** @todo Test Logarithmic interpolation: - * 1) how slow is it? - * 2) how much more accurate it is? (propably not much) - * 3) if not worth to do it replace with: - * F1 = F1F2DeviceBuff(DataSpace<2>{index, 0}); - * F2 = F1F2DeviceBuff(DataSpace<2>{index, 1}); - */ - - //! logInterpolate F1 and F2 -> passed by reference logInterpolation(index, F1, Ftemp, f, 0, F1F2DeviceBuff); logInterpolation(index, F2, Ftemp, f, 1, F1F2DeviceBuff); } - //! Calculating the numeric factor: dt * (e**2 * m_e * c /( hbar**2 * eps0 * 4 * np.pi)) - float_X numericFactor = DELTA_T + return std::make_tuple(F1, F2); + } + + //! Calculates the numeric factor for probability calculation + HDINLINE float_X calculateNumericFactor() + { + return DELTA_T * (ELECTRON_CHARGE * ELECTRON_CHARGE * ELECTRON_MASS * SPEED_OF_LIGHT / (HBAR * HBAR * EPS0 * 4._X * PI)); - + } + //! Checks if the probability calculation requirements are met and updates the buffer if not + template + HDINLINE void checkRequirements( + const T_Worker& worker, + float_X chi, + float_X gamma, + float_X numericFactor, + GridBuffer::DataBoxType failedRequirementQBuff) + { if constexpr(params::supressRequirementWarning == false) { float_X requirement1 = numericFactor * 1.5_X * math::pow(chi, 2._X / 3._X) / gamma; float_X requirement2 = numericFactor * 0.5_X * chi / gamma; - /** this warning means that the propability of generating a photon is high for given dt (higher - * than 10%) this means that we possibly generate photons every timestep and - * underestimate the radiation. The solution would be to implement substepping. - * If requirement is not met send warning at every iteration - */ if(requirement1 > 0.1_X || requirement2 > 0.1_X) + { alpaka::atomicExch(worker.getAcc(), &(failedRequirementQBuff(DataSpace<1>{0})), 1); + } } - - numericFactor *= math::sqrt(3._X) / (pmacc::math::Pi::doubleValue); - - //! Calculate propability: - float_X const numerator - = oneMinusDeltaOverDelta * chi * (F1 + 3._X * delta * zq * chi / 2._X * F2); - //! the delta is in oneMinusDeltaOverDelta in numerator1 - float_X const denominator = gamma; - - float_X propability = numericFactor * numerator / denominator * 3._X * r1r1; - - return (propability > randNr2) * delta; } - }; - - /** \struct AlgorithmSynchrotron - * This struct is passed to creation::createParticlesFromSpecies in the file: - * include/picongpu/particles/ParticlesFunctors.hpp - * - * @tparam T_DestSpecies type or name as PMACC_CSTRING of the photon species to be created - * @tparam T_SrcSpecies type or name as PMACC_CSTRING of the particle species that radiates so electrons - * only - * - * Takes care of: - * - random number generation - * - E and B field interpolation - * - maby something else as well - */ - template - struct AlgorithmSynchrotron - { - using DestSpecies = pmacc::particles::meta::FindByNameOrType_t; - using SrcSpecies = pmacc::particles::meta::FindByNameOrType_t; - - using FrameType = typename SrcSpecies::FrameType; - - /** specify field to particle interpolation scheme */ - using Field2ParticleInterpolation = typename pmacc::traits::Resolve< - typename pmacc::traits::GetFlagType>::type>::type; - - /** margins around the supercell for the interpolation of the field on the cells */ - using LowerMargin = typename GetMargin::LowerMargin; - using UpperMargin = typename GetMargin::UpperMargin; - - /** relevant area of a block */ - using BlockArea = SuperCellDescription; - - BlockArea BlockDescription; - - private: - /** random number generator */ - using RNGFactory = pmacc::random::RNGProvider; - using Distribution = pmacc::random::distributions::Uniform; - using RandomGen = typename RNGFactory::GetRandomType::type; - RandomGen randomGen; - - using TVec = MappingDesc::SuperCellSize; - - using ValueType_E = FieldE::ValueType; - using ValueType_B = FieldB::ValueType; - /** global memory EM-field and current density device databoxes */ - PMACC_ALIGN(eBox, FieldE::DataBoxType); - PMACC_ALIGN(bBox, FieldB::DataBoxType); - /** shared memory EM-field device databoxes */ - PMACC_ALIGN(cachedE, DataBox>); - PMACC_ALIGN(cachedB, DataBox>); - - //! F1F2DeviceBuff_ is a pointer to a databox containing F1 and F2 values. m stands for member - GridBuffer::DataBoxType m_F1F2DeviceBuff; - GridBuffer::DataBoxType m_failedRequirementQBuff; - - public: - /** host constructor initializing member : random number generator */ - AlgorithmSynchrotron( - const uint32_t currentStep, - GridBuffer::DataBoxType F1F2DeviceBuff, - GridBuffer::DataBoxType failedRequirementQBuff) - : randomGen(RNGFactory::createRandom()) - , m_F1F2DeviceBuff{F1F2DeviceBuff} - , m_failedRequirementQBuff{failedRequirementQBuff} + //! Calculates the probability of photon generation + HDINLINE float_X calculateProbability( + float_X numericFactor, + float_X oneMinusDeltaOverDelta, + float_X chi, + float_X zq, + float_X F1, + float_X F2, + float_X delta, + float_X gamma, + float_X r1r1, + float_64 randNr2) { - DataConnector& dc = Environment<>::get().DataConnector(); - /** initialize pointers on host-side E-(B-)field and current density databoxes */ - auto fieldE = dc.get(FieldE::getName()); - auto fieldB = dc.get(FieldB::getName()); - /** initialize device-side E-(B-)field and current density databoxes */ - eBox = fieldE->getDeviceDataBox(); - bBox = fieldB->getDeviceDataBox(); + numericFactor *= math::sqrt(3._X) / (pmacc::math::Pi::doubleValue); + float_X numerator = oneMinusDeltaOverDelta * chi * (F1 + 3._X * delta * zq * chi / 2._X * F2); + float_X denominator = gamma; + float_X probability = numericFactor * numerator / denominator * 3._X * r1r1; + return (probability > randNr2) * delta; } - /** cache fields used by this functor - * - * @warning this is a collective method and calls synchronize + /** Functor implementation + * @tparam EType type of electric field + * @tparam BType type of magnetic field + * @tparam ParticleType type of electron particle * - * @tparam T_Worker lockstep::Worker, lockstep worker + * @param bField magnetic field value at t=0 + * @param eField electric field value at t=0 + * @param parentElectron instance of electron with position at t=0 and momentum at t=-1/2 + * @param randNr random number, equally distributed in range [0.:1.0] * - * @param worker lockstep worker - * @param blockCell relative offset (in cells) to the local domain plus the guarding cells - * @param workerCfg configuration of the worker + * @return photon energy with respect to electron energy if photon is created 0 otherwise */ - template - DINLINE void collectiveInit(const T_Worker& worker, const DataSpace& blockCell) + template + HDINLINE float_X operator()( + const T_Worker& worker, + const BType bField, + const EType eField, + ParticleType& parentElectron, + float_64 randNr1, + float_64 randNr2, + GridBuffer::DataBoxType F1F2DeviceBuff, + GridBuffer::DataBoxType failedRequirementQBuff) const { - /** caching of E and B fields */ - cachedB = CachedBox::create<0, ValueType_B>(worker, BlockArea()); - cachedE = CachedBox::create<1, ValueType_E>(worker, BlockArea()); - - /** instance of alpaka assignment operator */ - pmacc::math::operation::Assign assign; - /** copy fields from global to shared */ - auto fieldBBlock = bBox.shift(blockCell); - auto collective = makeThreadCollective(); - collective(worker, assign, cachedB, fieldBBlock); - /** copy fields from global to shared */ - auto fieldEBlock = eBox.shift(blockCell); - collective(worker, assign, cachedE, fieldEBlock); - - /** wait for shared memory to be initialized */ - worker.sync(); - } + auto [minZqExponent, maxZqExponent, interpolationPoints, stepWidthLogarithmicScale] + = getInterpolationParams(); - /** Initialization function on device - * - * \brief Cache EM-fields on device - * and initialize possible prerequisites for radiation, 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 `cupla::__syncthreads( acc )` call which is - * necessary after initializing the E-/B-field shared boxes in shared memory. - */ - template - DINLINE void init( - T_Worker const&, - const DataSpace& localSuperCellOffset, - const uint32_t rngIdx) + float_X delta = calculateDelta(randNr1); + float_X HeffValue = calculateHeff(parentElectron, bField, eField); - { - auto rngOffset = DataSpace::create(0); - rngOffset.x() = rngIdx; - auto numRNGsPerSuperCell = DataSpace::create(1); - numRNGsPerSuperCell.x() = FrameType::frameSize; - this->randomGen.init(localSuperCellOffset * numRNGsPerSuperCell + rngOffset); - } + if(HeffValue == 0) + { + return 0; + } - /** Determine number of new macro photons due to radiation -> called by CreationKernel - * - * @param electronFrame reference to frame of the electron - * @param localIdx local (linear) index in super cell / frame - */ - template - DINLINE uint32_t numNewParticles(const T_Worker& worker, FrameType& electronFrame, int localIdx) - { - /** alias for the single macro-particle - electron */ - auto particle = electronFrame[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 = pmacc::math::mapToND(TVec::toRT(), particleCellIdx); - /** interpolation of E-... */ - const picongpu::traits::FieldPosition fieldPosE; - ValueType_E eField = Field2ParticleInterpolation()(cachedE.shift(localCell), pos, fieldPosE()); - /** ...and B-field on the particle position */ - const picongpu::traits::FieldPosition fieldPosB; - ValueType_B bField = Field2ParticleInterpolation()(cachedB.shift(localCell), pos, fieldPosB()); - - //! use the algorithm from the SynchrotronIdea struct - SynchrotronIdea synchrotronAlgo; - float_X photonEnergy = synchrotronAlgo( - worker, - bField, - eField, - particle, - this->randomGen(worker), - this->randomGen(worker), - m_F1F2DeviceBuff, - m_failedRequirementQBuff); - - - //! conversion factor from photon energy to momentum - constexpr float_X convFactor = 1.0_X / SPEED_OF_LIGHT; - //! if this is wrong uncomment the lines below and comment this line - float3_X const PhotonMomentum = particle[momentum_] * photonEnergy * convFactor; - - //! save to member variable to use in creation of new photon - m_PhotonMomentum = PhotonMomentum; - - //! generate one photon if energy is > minEnergy - return (pmacc::math::l2norm(PhotonMomentum) > params::minEnergy * convFactor); + float_X chi = calculateChi(parentElectron, HeffValue); + auto [zq, index] = calculateZq(chi, delta, minZqExponent, stepWidthLogarithmicScale); + auto [F1, F2] = interpolateF1F2( + index, + zq, + minZqExponent, + interpolationPoints, + stepWidthLogarithmicScale, + F1F2DeviceBuff); + + float_X numericFactor = calculateNumericFactor(); + checkRequirements(worker, chi, gamma, numericFactor, failedRequirementQBuff); + + return calculateProbability( + numericFactor, + oneMinusDeltaOverDelta, + chi, + zq, + F1, + F2, + delta, + gamma, + r1r1, + randNr2); } - /** Functor implementation + /** \struct AlgorithmSynchrotron + * This struct is passed to creation::createParticlesFromSpecies in the file: + * include/picongpu/particles/ParticlesFunctors.hpp * - * Ionization model specific particle creation + * @tparam T_DestSpecies type or name as PMACC_CSTRING of the photon species to be created + * @tparam T_SrcSpecies type or name as PMACC_CSTRING of the particle species that radiates so + * electrons only * - * @tparam T_parentElectron type of ion species that is radiating - * @tparam T_childPhoton type of Photon species that is created - * @param parentElectron electron instance that radiates - * @param childPhoton Photon instance that is created + * Takes care of: + * - random number generation + * - E and B field interpolation + * - maby something else as well */ - template - DINLINE void operator()( - const T_Worker& worker, - T_parentElectron& parentElectron, - T_childPhoton& childPhoton) + template + struct AlgorithmSynchrotron { - /** for not mixing operations::assign up with the nvidia functor assign */ - namespace partOp = pmacc::particles::operations; - /** each thread sets the multiMask hard on "particle" (=1) */ - childPhoton[multiMask_] = 1u; - - /** each thread initializes a clone of the parent Electron but leaving out - * some attributes: - * - multiMask: reading from global memory takes longer than just setting it again explicitly - * - momentum: because the Photon would get a higher energy because of the Electron mass + using DestSpecies = pmacc::particles::meta::FindByNameOrType_t; + using SrcSpecies = pmacc::particles::meta::FindByNameOrType_t; + + using FrameType = typename SrcSpecies::FrameType; + + /** specify field to particle interpolation scheme */ + using Field2ParticleInterpolation = typename pmacc::traits::Resolve< + typename pmacc::traits::GetFlagType>::type>::type; + + /** margins around the supercell for the interpolation of the field on the cells */ + using LowerMargin = typename GetMargin::LowerMargin; + using UpperMargin = typename GetMargin::UpperMargin; + + /** relevant area of a block */ + using BlockArea + = SuperCellDescription; + + BlockArea BlockDescription; + + private: + /** random number generator */ + using RNGFactory = pmacc::random::RNGProvider; + using Distribution = pmacc::random::distributions::Uniform; + using RandomGen = typename RNGFactory::GetRandomType::type; + RandomGen randomGen; + + using TVec = MappingDesc::SuperCellSize; + + using ValueType_E = FieldE::ValueType; + using ValueType_B = FieldB::ValueType; + /** global memory EM-field and current density device databoxes */ + PMACC_ALIGN(eBox, FieldE::DataBoxType); + PMACC_ALIGN(bBox, FieldB::DataBoxType); + /** shared memory EM-field device databoxes */ + PMACC_ALIGN(cachedE, DataBox>); + PMACC_ALIGN(cachedB, DataBox>); + + //! F1F2DeviceBuff_ is a pointer to a databox containing F1 and F2 values. m stands for member + GridBuffer::DataBoxType m_F1F2DeviceBuff; + GridBuffer::DataBoxType m_failedRequirementQBuff; + + public: + /** host constructor initializing member : random number generator */ + AlgorithmSynchrotron( + const uint32_t currentStep, + GridBuffer::DataBoxType F1F2DeviceBuff, + GridBuffer::DataBoxType failedRequirementQBuff) + : randomGen(RNGFactory::createRandom()) + , m_F1F2DeviceBuff{F1F2DeviceBuff} + , m_failedRequirementQBuff{failedRequirementQBuff} + { + DataConnector& dc = Environment<>::get().DataConnector(); + /** initialize pointers on host-side E-(B-)field and current density databoxes */ + auto fieldE = dc.get(FieldE::getName()); + auto fieldB = dc.get(FieldB::getName()); + /** initialize device-side E-(B-)field and current density databoxes */ + eBox = fieldE->getDeviceDataBox(); + bBox = fieldB->getDeviceDataBox(); + } + + /** cache fields used by this functor + * + * @warning this is a collective method and calls synchronize + * + * @tparam T_Worker lockstep::Worker, lockstep worker + * + * @param worker lockstep worker + * @param blockCell relative offset (in cells) to the local domain plus the guarding cells + * @param workerCfg configuration of the worker + */ + template + DINLINE void collectiveInit(const T_Worker& worker, const DataSpace& blockCell) + { + /** caching of E and B fields */ + cachedB = CachedBox::create<0, ValueType_B>(worker, BlockArea()); + cachedE = CachedBox::create<1, ValueType_E>(worker, BlockArea()); + + /** instance of alpaka assignment operator */ + pmacc::math::operation::Assign assign; + /** copy fields from global to shared */ + auto fieldBBlock = bBox.shift(blockCell); + auto collective = makeThreadCollective(); + collective(worker, assign, cachedB, fieldBBlock); + /** copy fields from global to shared */ + auto fieldEBlock = eBox.shift(blockCell); + collective(worker, assign, cachedE, fieldEBlock); + + /** wait for shared memory to be initialized */ + worker.sync(); + } + + /** Initialization function on device + * + * \brief Cache EM-fields on device + * and initialize possible prerequisites for radiation, 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 `cupla::__syncthreads( acc )` call which is + * necessary after initializing the E-/B-field shared boxes in shared memory. */ - auto targetPhotonClone = partOp::deselect>(childPhoton); + template + DINLINE void init( + T_Worker const&, + const DataSpace& localSuperCellOffset, + const uint32_t rngIdx) - partOp::assign(targetPhotonClone, partOp::deselect(parentElectron)); + { + auto rngOffset = DataSpace::create(0); + rngOffset.x() = rngIdx; + auto numRNGsPerSuperCell = DataSpace::create(1); + numRNGsPerSuperCell.x() = FrameType::frameSize; + this->randomGen.init(localSuperCellOffset * numRNGsPerSuperCell + rngOffset); + } - childPhoton[momentum_] = m_PhotonMomentum; + /** Determine number of new macro photons due to radiation -> called by CreationKernel + * + * @param electronFrame reference to frame of the electron + * @param localIdx local (linear) index in super cell / frame + */ + template + DINLINE uint32_t numNewParticles(const T_Worker& worker, FrameType& electronFrame, int localIdx) + { + /** alias for the single macro-particle - electron */ + auto particle = electronFrame[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 = pmacc::math::mapToND(TVec::toRT(), particleCellIdx); + /** interpolation of E-... */ + const picongpu::traits::FieldPosition fieldPosE; + ValueType_E eField = Field2ParticleInterpolation()(cachedE.shift(localCell), pos, fieldPosE()); + /** ...and B-field on the particle position */ + const picongpu::traits::FieldPosition fieldPosB; + ValueType_B bField = Field2ParticleInterpolation()(cachedB.shift(localCell), pos, fieldPosB()); + + //! use the algorithm from the SynchrotronIdea struct + SynchrotronIdea synchrotronAlgo; + float_X photonEnergy = synchrotronAlgo( + worker, + bField, + eField, + particle, + this->randomGen(worker), + this->randomGen(worker), + m_F1F2DeviceBuff, + m_failedRequirementQBuff); + + + //! conversion factor from photon energy to momentum + constexpr float_X convFactor = 1.0_X / SPEED_OF_LIGHT; + //! if this is wrong uncomment the lines below and comment this line + float3_X const PhotonMomentum = particle[momentum_] * photonEnergy * convFactor; + + //! save to member variable to use in creation of new photon + m_PhotonMomentum = PhotonMomentum; + + //! generate one photon if energy is > minEnergy + return (pmacc::math::l2norm(PhotonMomentum) > params::minEnergy * convFactor); + } - /** conservatElectron of momentum */ - if constexpr(params::ElectronRecoil) - parentElectron[momentum_] -= m_PhotonMomentum; - } + /** Functor implementation + * + * Ionization model specific particle creation + * + * @tparam T_parentElectron type of ion species that is radiating + * @tparam T_childPhoton type of Photon species that is created + * @param parentElectron electron instance that radiates + * @param childPhoton Photon instance that is created + */ + template + DINLINE void operator()( + const T_Worker& worker, + T_parentElectron& parentElectron, + T_childPhoton& childPhoton) + { + /** for not mixing operations::assign up with the nvidia functor assign */ + namespace partOp = pmacc::particles::operations; + /** each thread sets the multiMask hard on "particle" (=1) */ + childPhoton[multiMask_] = 1u; + + /** each thread initializes a clone of the parent Electron but leaving out + * some attributes: + * - multiMask: reading from global memory takes longer than just setting it again explicitly + * - momentum: because the Photon would get a higher energy because of the Electron mass + */ + auto targetPhotonClone = partOp::deselect>(childPhoton); + + partOp::assign(targetPhotonClone, partOp::deselect(parentElectron)); + + childPhoton[momentum_] = m_PhotonMomentum; + + /** conservatElectron of momentum */ + if constexpr(params::ElectronRecoil) + parentElectron[momentum_] -= m_PhotonMomentum; + } - private: - //! energy of emitted photon with respect to electron energy - float3_X m_PhotonMomentum; + private: + //! energy of emitted photon with respect to electron energy + float3_X m_PhotonMomentum; - }; // end of struct AlgorithmSynchrotron - } // namespace synchrotron - } // namespace particles -} // namespace picongpu + }; // end of struct AlgorithmSynchrotron + } // namespace synchrotron + } // namespace particles + } // namespace picongpu From 432b700610d7446934e418d13754b90be91bf642 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Filip=20Opto=C5=82owicz?= Date: Fri, 7 Jun 2024 17:51:21 +0200 Subject: [PATCH 2/2] some bug fixes ideas --- .../synchrotron/AlgorithmSynchrotron.hpp | 26 ++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp b/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp index 054297fe9a..e042610635 100644 --- a/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp +++ b/include/picongpu/particles/synchrotron/AlgorithmSynchrotron.hpp @@ -437,7 +437,31 @@ namespace picongpu //! conversion factor from photon energy to momentum constexpr float_X convFactor = 1.0_X / SPEED_OF_LIGHT; //! if this is wrong uncomment the lines below and comment this line - float3_X const PhotonMomentum = particle[momentum_] * photonEnergy * convFactor; + float3_X const ratioPhotonElectronEnergy = particle[momentum_] * photonEnergy * convFactor; + + //! conversion factor from photon energy to momentum + constexpr float_X convFactor = 1.0_X / SPEED_OF_LIGHT; + // E*r = E_y = p_y*c + // + // E^2 = p^2 * c^2 + m^2 * c^4 + // p^2 = (E^2 - m^2 * c^4) / c^2 + // + // p_y*c = sqrt((p^2 * c^2 + m^2 * c^4)*r^2) + // p_y = sqrt((p^2 + m^2 * c^2) * r^2) + // p_y = sqrt((p^2 * r^2 + m^2 * c^2 * r^2)) != p*r/c + + //! if this is wrong uncomment the lines below and comment this line + // float3_X const PhotonMomentum = particle[momentum_] * ratioPhotonElectronMomenta; + + float3_X const mom = particle[momentum_]; + float_X const mass = attribute::getMass(float_X(1), particle); // weighting 1 + + float_X const gamma = Gamma()(mom, mass); + constexpr float_X c2 = SPEED_OF_LIGHT * SPEED_OF_LIGHT; + float_X const electronKineticEnergy = (gamma - float_X(1.0)) * mass * c2; + float_X const photonEnergy = electronKineticEnergy * ratioPhotonElectronEnergy; + float3_X const PhotonMomentum = mom / pmacc::math::l2norm(mom) * photonEnergy * convFactor; + //! save to member variable to use in creation of new photon m_PhotonMomentum = PhotonMomentum;