diff --git a/include/picongpu/particles/manipulators/manipulators.def b/include/picongpu/particles/manipulators/manipulators.def index 431989fb7f..d6861e1b2c 100644 --- a/include/picongpu/particles/manipulators/manipulators.def +++ b/include/picongpu/particles/manipulators/manipulators.def @@ -35,6 +35,7 @@ #include "picongpu/particles/manipulators/unary/Drift.def" #include "picongpu/particles/manipulators/unary/FreeTotalCellOffset.def" #include "picongpu/particles/manipulators/unary/FreeTotalCellOffsetRng.def" +#include "picongpu/particles/manipulators/unary/MaxwellJuettner.def" #include "picongpu/particles/manipulators/unary/None.def" #include "picongpu/particles/manipulators/unary/RandomPosition.def" #include "picongpu/particles/manipulators/unary/SparserMacroParticlesPerSuperCell.def" diff --git a/include/picongpu/particles/manipulators/manipulators.hpp b/include/picongpu/particles/manipulators/manipulators.hpp index 2bdd7a4bbe..065f3b9c7c 100644 --- a/include/picongpu/particles/manipulators/manipulators.hpp +++ b/include/picongpu/particles/manipulators/manipulators.hpp @@ -25,5 +25,6 @@ #include "picongpu/particles/manipulators/unary/Drift.hpp" #include "picongpu/particles/manipulators/unary/FreeTotalCellOffset.hpp" #include "picongpu/particles/manipulators/unary/FreeTotalCellOffsetRng.hpp" +#include "picongpu/particles/manipulators/unary/MaxwellJuettner.hpp" #include "picongpu/particles/manipulators/unary/SparserMacroParticlesPerSuperCell.hpp" #include "picongpu/particles/manipulators/unary/Temperature.hpp" diff --git a/include/picongpu/particles/manipulators/unary/MaxwellJuettner.def b/include/picongpu/particles/manipulators/unary/MaxwellJuettner.def new file mode 100644 index 0000000000..a404ba4ade --- /dev/null +++ b/include/picongpu/particles/manipulators/unary/MaxwellJuettner.def @@ -0,0 +1,156 @@ +/* Copyright 2014-2024 Rene Widera, Sergei Bastrakov, Sergey Ermakov, + * Sergey Ermakov + * + * 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/simulation_defines.hpp" + +#include "picongpu/particles/manipulators/generic/FreeRng.def" +#include "picongpu/particles/manipulators/unary/FreeTotalCellOffsetRng.def" + +#include +#include + +#include + + +namespace picongpu::particles::manipulators::unary +{ + namespace acc + { + /** Functor to modify particle momentum based on temperature for relativistic case + * Use Maxwell Juettner Distribution, implementation based on https://doi.org/10.1063/1.4919383 + * This functor is for the highly-relativistic case. + * In this case the added momentum follows the Maxwell-Juettner distribution. + * + * @tparam T_ValueFunctor pmacc::math::operation::*, binary functor type to + * add a new momentum to an old one + */ + + /** Version with fixed temperature given via parameter struct + * + * @tparam T_ParamClass configuration parameter, follows requirements of param::MaxwellJuettnerCfg + * + * USE THIS IN THE SAME WAY AS Temperature.hpp JUST REPLACE Temperature with MaxwellJuettner FOR + * ALL FUNCTOR AND TEMPLATE NAMES + **/ + template + struct MaxwellJuettner; + + /** Version with user-provided temperature functor + * + * @tparam T_MaxwellJuettnerFunctor temperature functor, follows requirements of + * param::MaxwellJuettnerFunctor + **/ + template + struct FreeMaxwellJuettner; + + /** @} */ + } // namespace acc + + namespace param + { + //! Configuration struct for the unary manipulator MaxwellJuettner + struct MaxwellJuettnerCfg + { + /** Initial temperature + * unit: keV + */ + static constexpr float_64 temperature = 0.0; + }; + + //! Configuration functor for the unary manipulator FreeMaxwellJuettner + struct MaxwellJuettnerFunctor + { + /** Constructors + * + * At least one of them must be provided (or auto-generated by a compiler). + * The version with current step will be called when available. + * + * @{ + */ + + //! Default constructor + HINLINE MaxwellJuettnerFunctor() + { + } + + /** Constructor with current step + * + * @param currentStep current time iteration + */ + HINLINE MaxwellJuettnerFunctor(uint32_t currentStep) + { + } + + /** @} */ + + /** Return the temperature in keV for the given position + * + * Return type may be float_X or float_64. + * + * @param position_SI total offset including all slides [meter] + * @param cellSize_SI cell sizes [meter] + */ + HDINLINE float_X operator()(floatD_64 const& position_SI, float3_64 const& cellSize_SI) + { + return 0.0_X; + } + }; + } // namespace param + + /** Modify particle momentum based on temperature + * + * Sample a random momentum value distributed according to the given temperature and add it to the + * existing particle momentum. + * + * Note: initial electron temperature should generally be chosen so that the resulting Debye length is + * resolved by the grid. + * + * @tparam T_ValueFunctor pmacc::math::operation::*, binary functor type to + * add a new momentum to an old one + * + * @{ + */ + + /** Version with fixed temperature given via parameter struct + * + * @tparam T_ParamClass configuration parameter, follows requirements of param::MaxwellJuettnerCfg + */ + template + using MaxwellJuettner = generic::FreeRng< + acc::MaxwellJuettner, + pmacc::random::distributions::Uniform>>; + + /** Version with user-provided temperature functor + * + * @tparam T_MaxwellJuettnerFunctor temperature functor, follows requirements of + * param::MaxwellJuettnerFunctor + */ + template< + typename T_MaxwellJuettnerFunctor = param::MaxwellJuettnerFunctor, + typename T_ValueFunctor = pmacc::math::operation::Add> + using FreeMaxwellJuettner = FreeTotalCellOffsetRng< + acc::FreeMaxwellJuettner, + pmacc::random::distributions::Uniform>>; + + /** @} */ + +} // namespace picongpu::particles::manipulators::unary diff --git a/include/picongpu/particles/manipulators/unary/MaxwellJuettner.hpp b/include/picongpu/particles/manipulators/unary/MaxwellJuettner.hpp new file mode 100644 index 0000000000..5c80281f7c --- /dev/null +++ b/include/picongpu/particles/manipulators/unary/MaxwellJuettner.hpp @@ -0,0 +1,176 @@ +/* Copyright 2013-2024 Axel Huebl, Heiko Burau, Rene Widera, + * Alexander Grund, Sergei Bastrakov, Sergey Ermakov + * + * 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/simulation_defines.hpp" + +#include "picongpu/particles/functor/User.hpp" + +namespace picongpu::particles::manipulators::unary::acc +{ + namespace detail + { + /** Functor to modify particle momentum based on temperature for relativistic case + * Use Maxwell Juettner Distribution, implementation based on https://doi.org/10.1063/1.4919383 + * This functor is for the highly-relativistic case. + * In this case the added momentum follows the Maxwell-Juettner distribution. + * + * @tparam T_ValueFunctor pmacc::math::operation::*, binary functor type to + * add a new momentum to an old one + */ + /** USE THIS IN THE SAME WAY AS Temperature.hpp JUST REPLACE Temperature with MaxwellJuettner + * FOR ALL FUNCTOR AND TEMPLATE NAMES**/ + template + struct MaxwellJuettnerImpl : private T_ValueFunctor + { + /** Manipulate the momentum of the given macroparticle + * + * @tparam T_UniformRng functor::misc::RngWrapper, standard + * normal random number generator type + * @tparam T_Particle particle type + * @tparam T_MaxwellJuettner temperature type + * + * @param uniformRng random number generator + * @param particle particle to be manipulated + * @param temperatureKeV temperature value in keV + */ + template + HDINLINE void operator()( + T_UniformNoZeroRng& uniformNoZeroRng, + T_Particle& particle, + T_MaxwellJuettner temperatureKeV) const + { + float_X const energy = (temperatureKeV * sim.si.conv.ev2Joule(1.0e3)) / sim.unit.energy(); + + float_X const macroWeighting = particle[weighting_]; + float_X const macroMass = attribute::getMass(macroWeighting, particle); + float_X const mass = macroMass / macroWeighting; + float_X const theta = energy / (mass * sim.pic.getSpeedOfLight() * sim.pic.getSpeedOfLight()); + + float_X x1, x2, x3, x4, u, nu; + bool rejected = true; + while(rejected) + { + x1 = uniformNoZeroRng(); + x2 = uniformNoZeroRng(); + x3 = uniformNoZeroRng(); + + u = (-1._X) * theta * math::log(x1 * x2 * x3); + + x4 = uniformNoZeroRng(); + nu = (-1._X) * theta * math::log(x1 * x2 * x3 * x4); + + if(nu * nu - u * u > 1._X) + { + rejected = false; + } + } + + float_X momAbs = u * macroWeighting * mass; + float_X y1 = uniformNoZeroRng(); + float_X y2 = uniformNoZeroRng(); + + float_X sin_2_pi_y2; + float_X cos_2_pi_y2; + math::sincos(precisionCast(2._X * PI * y2), sin_2_pi_y2, cos_2_pi_y2); + + auto mom = float3_X( + momAbs * (2._X * y1 - 1._X), + 2._X * momAbs * math::sqrt(y1 * (1._X - y1)) * cos_2_pi_y2, + 2._X * momAbs * math::sqrt(y1 * (1._X - y1)) * sin_2_pi_y2); + + T_ValueFunctor::operator()(particle[momentum_], mom); + } + }; + + } // namespace detail + + template + struct MaxwellJuettner : public detail::MaxwellJuettnerImpl + { + //! Base class + using Base = detail::MaxwellJuettnerImpl; + + /** Manipulate the momentum of the given macroparticle + * + * @tparam T_UniformRng functor::misc::RngWrapper, + * uniform random number generator type + * @tparam T_Particle particle type + * + * @param uniform random number generator + * @param particle particle to be manipulated + */ + template + HDINLINE void operator()(T_UniformRng& uniformRng, T_Particle& particle) + { + auto const temperatureKeV = T_ParamClass::temperature; + Base::operator()(uniformRng, particle, temperatureKeV); + } + }; + + template + struct FreeMaxwellJuettner + : public detail::MaxwellJuettnerImpl + , public particles::functor::User + { + //! Base implementation class + using Base = detail::MaxwellJuettnerImpl; + + //! Wrapper around user-provided functor + using UserFunctor = particles::functor::User; + + /** Create a functor instance, including instances for user functor and its wrapper + * + * @param currentStep current time iteration + */ + FreeMaxwellJuettner(uint32_t const currentStep, IdGenerator idGen) + : UserFunctor(currentStep, idGen) + { + } + + /** Manipulate the momentum of the given macroparticle + * + * @tparam T_UniformRng functor::misc::RngWrapper, standard + * uniform random number generator type + * @tparam T_Particle particle type + * @tparam T_Args arbitrary number of argument types, unused + * + * @param totalCellOffset total offset including all slides [in cells] + * @param uniformRng uniform random number generator + * @param particle particle to be manipulated + * @param ... unused parameters + */ + template + HDINLINE void operator()( + DataSpace const& totalCellOffset, + T_UniformRng& uniformRng, + T_Particle& particle) + { + auto const unitLength = sim.unit.length(); + auto const cellSize_SI = precisionCast(sim.pic.getCellSize()) * unitLength; + auto const position_SI + = (precisionCast(totalCellOffset) + precisionCast(particle[position_])) + * cellSize_SI.shrink(); + auto const temperatureKeV = UserFunctor::operator()(position_SI, cellSize_SI); + Base::operator()(uniformRng, particle, temperatureKeV); + } + }; +} // namespace picongpu::particles::manipulators::unary::acc