diff --git a/include/picongpu/param/pusher.param b/include/picongpu/param/pusher.param index dcf5597be0..a30fa58e93 100644 --- a/include/picongpu/param/pusher.param +++ b/include/picongpu/param/pusher.param @@ -62,6 +62,15 @@ namespace picongpu namespace sqrt_Vay = precision64Bit; } + namespace particlePusherHigueraCary + { + /** Precision of the square roots during the push step + * - precision32Bit + * - precision64Bit + */ + namespace sqrt_HigueraCary = precision64Bit; + } + namespace particlePusherAxel { @@ -78,6 +87,7 @@ namespace picongpu { namespace pusher { + struct HigueraCary; struct Vay; struct Boris; struct Photon; diff --git a/include/picongpu/param/species.param b/include/picongpu/param/species.param index 5ae23aa732..5d1a6b3f06 100644 --- a/include/picongpu/param/species.param +++ b/include/picongpu/param/species.param @@ -1,4 +1,4 @@ -/* Copyright 2014-2020 Rene Widera, Richard Pausch +/* Copyright 2014-2020 Rene Widera, Richard Pausch, Annegret Roeszler * * This file is part of PIConGPU. * @@ -74,8 +74,9 @@ using UsedParticleCurrentSolver = currentSolver::Esirkepov< UsedParticleShape >; * * Defining a pusher is optional for particles * - * - particles::pusher::Vay : better suited relativistic boris pusher - * - particles::pusher::Boris : standard boris pusher + * - particles::pusher::HigueraCary : Higuera & Cary's relativistic pusher preserving both volume and ExB velocity + * - particles::pusher::Vay : Vay's relativistic pusher preserving ExB velocity + * - particles::pusher::Boris : Boris' relativistic pusher preserving volume * - particles::pusher::ReducedLandauLifshitz : 4th order RungeKutta pusher * with classical radiation reaction * diff --git a/include/picongpu/particles/pusher/particlePusherHigueraCary.hpp b/include/picongpu/particles/pusher/particlePusherHigueraCary.hpp new file mode 100644 index 0000000000..f33c924fe3 --- /dev/null +++ b/include/picongpu/particles/pusher/particlePusherHigueraCary.hpp @@ -0,0 +1,145 @@ +/* Copyright 2013-2020 Heiko Burau, Rene Widera, Richard Pausch, Annegret Roeszler, Klaus Steiniger + * + * 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/traits/attribute/GetMass.hpp" +#include "picongpu/traits/attribute/GetCharge.hpp" + + +namespace picongpu +{ +namespace particlePusherHigueraCary +{ + +/** Implementation of the Higuera-Cary pusher as presented in doi:10.1063/1.4979989. + * + * A correction is applied to the given formulas as documented by the WarpX team: + * (https://github.com/ECP-WarpX/WarpX/issues/320). + * + * Note, while Higuera and Ripperda present the formulas for the quantity u = gamma * v, + * PIConGPU uses the real momentum p = gamma * m * v = u * m for calculations. + * Here, all auxiliary quantities are equal to those used in Ripperda's article. + * + * Further references: + * [Higuera's article on arxiv](https://arxiv.org/abs/1701.05605) + * [Riperda's comparison of relativistic particle integrators](https://doi.org/10.3847/1538-4365/aab114) + * + * @tparam Velocity functor to compute the velocity of a particle with momentum p and mass m + * @tparam Gamma functor to compute the Lorentz factor (= Energy/mc^2) of a particle with momentum p and mass m + */ +template< + typename Velocity, + typename Gamma +> +struct Push +{ + /* this is an optional extension for sub-sampling pushes that enables grid to particle interpolation + * for particle positions outside the super cell in one push + */ + using LowerMargin = typename pmacc::math::CT::make_Int::type; + using UpperMargin = typename pmacc::math::CT::make_Int::type; + + template< typename T_FunctorFieldE, typename T_FunctorFieldB, typename T_Particle, typename T_Pos > + HDINLINE void operator()( + T_FunctorFieldB const functorBField, + T_FunctorFieldE const functorEField, + T_Particle & particle, + T_Pos & pos, + uint32_t const + ) + { + float_X const weighting = particle[ weighting_ ]; + float_X const mass = attribute::getMass( weighting , particle ); + float_X const charge = attribute::getCharge( weighting , particle ); + + using MomType = momentum::type; + MomType const mom = particle[ momentum_ ]; + + auto bField = functorBField(pos); + auto eField = functorEField(pos); + + float_X const deltaT = DELTA_T; + + + Gamma gamma; + + /* Momentum update + * Notation is according to Ripperda's paper + */ + // First half electric field acceleration + namespace sqrt_HC = sqrt_HigueraCary; + + sqrt_HC::float3_X const mom_minus = precisionCast( mom + float_X(0.5) * charge * eField * deltaT ); + + // Auxiliary quantitites + sqrt_HC::float_X const gamma_minus = gamma( mom_minus , mass ); + + sqrt_HC::float3_X const tau = precisionCast( float_X(0.5) * bField * charge * deltaT / mass ); + + sqrt_HC::float_X const sigma = pmacc::math::abs2( gamma_minus ) - pmacc::math::abs2( tau ); + + sqrt_HC::float_X const u_star = pmacc::math::dot( mom_minus, tau ) / precisionCast( mass * SPEED_OF_LIGHT ); + + sqrt_HC::float_X const gamma_plus = math::sqrt( + sqrt_HC::float_X(0.5) * ( sigma + math::sqrt( + pmacc::math::abs2( sigma ) + sqrt_HC::float_X(4.0) * ( pmacc::math::abs2( tau ) + pmacc::math::abs2( u_star ) ) + ) ) + ); + + sqrt_HC::float3_X const t_vector = tau / gamma_plus; + + sqrt_HC::float_X const s = sqrt_HC::float_X(1.0) / ( sqrt_HC::float_X(1.0) + pmacc::math::abs2( t_vector ) ); + + // Rotation step + sqrt_HC::float3_X const mom_plus = s * ( mom_minus + + pmacc::math::dot( mom_minus , t_vector ) * t_vector + + pmacc::math::cross( mom_minus , t_vector ) + ); + + // Second half electric field acceleration (Note correction mom_minus -> mom_plus here compared to Ripperda) + MomType const mom_diff1 = float_X(0.5) * charge * eField * deltaT; + MomType const mom_diff2 = precisionCast( pmacc::math::cross( mom_plus , t_vector ) ); + MomType const mom_diff = mom_diff1 + mom_diff2; + + MomType const new_mom = precisionCast( mom_plus ) + mom_diff; + + particle[ momentum_ ] = new_mom; + + // Position update + Velocity velocity; + + float3_X const vel = velocity( new_mom , mass ); + + for( uint32_t d=0 ; d > { }; +struct HigueraCary : +public particlePusherHigueraCary::Push > +{ +}; + struct Free : public particlePusherFree::Push > { diff --git a/share/picongpu/examples/Bunch/include/picongpu/param/species.param b/share/picongpu/examples/Bunch/include/picongpu/param/species.param index 86aea1bffb..4f6e2fb152 100644 --- a/share/picongpu/examples/Bunch/include/picongpu/param/species.param +++ b/share/picongpu/examples/Bunch/include/picongpu/param/species.param @@ -65,8 +65,9 @@ using UsedParticleCurrentSolver = currentSolver::Esirkepov< UsedParticleShape >; * * Defining a pusher is optional for particles * - * - particles::pusher::Vay : better suited relativistic boris pusher - * - particles::pusher::Boris : standard boris pusher + * - particles::pusher::HigueraCary : Higuera & Cary's relativistic pusher preserving both volume and ExB velocity + * - particles::pusher::Vay : Vay's relativistic pusher preserving ExB velocity + * - particles::pusher::Boris : Boris' relativistic pusher preserving volume * - particles::pusher::ReducedLandauLifshitz : 4th order RungeKutta pusher * with classical radiation reaction * diff --git a/share/picongpu/examples/KelvinHelmholtz/include/picongpu/param/species.param b/share/picongpu/examples/KelvinHelmholtz/include/picongpu/param/species.param index 388415199c..fa2dedad45 100644 --- a/share/picongpu/examples/KelvinHelmholtz/include/picongpu/param/species.param +++ b/share/picongpu/examples/KelvinHelmholtz/include/picongpu/param/species.param @@ -71,8 +71,9 @@ using UsedParticleCurrentSolver = currentSolver::PARAM_CURRENTSOLVER; * * Defining a pusher is optional for particles * - * - particles::pusher::Vay : better suited relativistic boris pusher - * - particles::pusher::Boris : standard boris pusher + * - particles::pusher::HigueraCary : Higuera & Cary's relativistic pusher preserving both volume and ExB velocity + * - particles::pusher::Vay : Vay's relativistic pusher preserving ExB velocity + * - particles::pusher::Boris : Boris' relativistic pusher preserving volume * - particles::pusher::ReducedLandauLifshitz : 4th order RungeKutta pusher * with classical radiation reaction * diff --git a/share/picongpu/examples/LaserWakefield/include/picongpu/param/species.param b/share/picongpu/examples/LaserWakefield/include/picongpu/param/species.param index 6def4fdd61..398bd46adc 100644 --- a/share/picongpu/examples/LaserWakefield/include/picongpu/param/species.param +++ b/share/picongpu/examples/LaserWakefield/include/picongpu/param/species.param @@ -71,8 +71,9 @@ using UsedParticleCurrentSolver = currentSolver::PARAM_CURRENTSOLVER< UsedPartic * * Defining a pusher is optional for particles * - * - particles::pusher::Vay : better suited relativistic boris pusher - * - particles::pusher::Boris : standard boris pusher + * - particles::pusher::HigueraCary : Higuera & Cary's relativistic pusher preserving both volume and ExB velocity + * - particles::pusher::Vay : Vay's relativistic pusher preserving ExB velocity + * - particles::pusher::Boris : Boris' relativistic pusher preserving volume * - particles::pusher::ReducedLandauLifshitz : 4th order RungeKutta pusher * with classical radiation reaction * diff --git a/share/picongpu/examples/SingleParticleTest/include/picongpu/param/species.param b/share/picongpu/examples/SingleParticleTest/include/picongpu/param/species.param index d6765d8bc3..cf073af2c1 100644 --- a/share/picongpu/examples/SingleParticleTest/include/picongpu/param/species.param +++ b/share/picongpu/examples/SingleParticleTest/include/picongpu/param/species.param @@ -65,8 +65,9 @@ using UsedParticleCurrentSolver = currentSolver::Esirkepov; * * Defining a pusher is optional for particles * - * - particles::pusher::Vay : better suited relativistic boris pusher - * - particles::pusher::Boris : standard boris pusher + * - particles::pusher::HigueraCary : Higuera & Cary's relativistic pusher preserving both volume and ExB velocity + * - particles::pusher::Vay : Vay's relativistic pusher preserving ExB velocity + * - particles::pusher::Boris : Boris' relativistic pusher preserving volume * - particles::pusher::ReducedLandauLifshitz : 4th order RungeKutta pusher * with classical radiation reaction * diff --git a/share/picongpu/examples/TransitionRadiation/include/picongpu/param/species.param b/share/picongpu/examples/TransitionRadiation/include/picongpu/param/species.param index 86aea1bffb..4f6e2fb152 100644 --- a/share/picongpu/examples/TransitionRadiation/include/picongpu/param/species.param +++ b/share/picongpu/examples/TransitionRadiation/include/picongpu/param/species.param @@ -65,8 +65,9 @@ using UsedParticleCurrentSolver = currentSolver::Esirkepov< UsedParticleShape >; * * Defining a pusher is optional for particles * - * - particles::pusher::Vay : better suited relativistic boris pusher - * - particles::pusher::Boris : standard boris pusher + * - particles::pusher::HigueraCary : Higuera & Cary's relativistic pusher preserving both volume and ExB velocity + * - particles::pusher::Vay : Vay's relativistic pusher preserving ExB velocity + * - particles::pusher::Boris : Boris' relativistic pusher preserving volume * - particles::pusher::ReducedLandauLifshitz : 4th order RungeKutta pusher * with classical radiation reaction *