aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Docs/source/running_cpp/parameters.rst6
-rwxr-xr-xExamples/Tests/radiation_reaction/test_const_B_analytical/check.py188
-rw-r--r--Examples/Tests/radiation_reaction/test_const_B_analytical/inputs67
-rw-r--r--Regression/WarpX-tests.ini15
-rw-r--r--Source/Particles/PhotonParticleContainer.H6
-rw-r--r--Source/Particles/PhysicalParticleContainer.H8
-rw-r--r--Source/Particles/PhysicalParticleContainer.cpp84
-rw-r--r--Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H84
-rw-r--r--Source/Particles/RigidInjectedParticleContainer.cpp30
9 files changed, 477 insertions, 11 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst
index 09d639728..e77ca553e 100644
--- a/Docs/source/running_cpp/parameters.rst
+++ b/Docs/source/running_cpp/parameters.rst
@@ -344,6 +344,7 @@ Particle initialization
* ``ux`` ``uy`` ``uz`` for the particle momentum,
* ``Ex`` ``Ey`` ``Ez`` for the electric field on particles,
* ``Bx`` ``By`` ``Bz`` for the magnetic field on particles.
+
The particle positions are always included. Use
``<species>.plot_vars = none`` to plot no particle data, except
particle position.
@@ -373,6 +374,11 @@ Particle initialization
species (must be smaller than the atomic number of chemical element given
in `physical_element`).
+* ``<species>.do_classical_radiation_reaction`` (`int`) optional (default `0`)
+ Enables Radiation Reaction (or Radiation Friction) for the species. Species
+ must be either electrons or positrons. Boris pusher must be used for the
+ simulation
+
* ``<species>.do_qed`` (`int`) optional (default `0`)
If `<species>.do_qed = 0` all the QED effects are disabled for this species.
If `<species>.do_qed = 1` QED effects can be enabled for this species (see below)
diff --git a/Examples/Tests/radiation_reaction/test_const_B_analytical/check.py b/Examples/Tests/radiation_reaction/test_const_B_analytical/check.py
new file mode 100755
index 000000000..6d0ca4502
--- /dev/null
+++ b/Examples/Tests/radiation_reaction/test_const_B_analytical/check.py
@@ -0,0 +1,188 @@
+#! /usr/bin/env python3
+
+# This script contains few simple tests for the radiation reaction pusher
+# It initializes an electron or a positron with normalized momentum in different
+# directions, propagating in a static magnetic field (along [2/7,3/7,6/7]).
+# If the initial momentum is perpendicular to the the field, there is a very
+# simple analytical formula for the gamma factor:
+# gamma(t) = coth(t/tc - C)
+# where tc = 1./(t0 * omega_c * omega_c)
+# and omega_c = qB/m , t0 = (2/3) * (r_e / c)
+# and r_e = classical electron radius
+# C is chosen so that gamma(0) = initial_gamma
+# If the initial momentum is parallel to the field, it should not change.
+# This test checks these predictions with a tolerance of 5%.
+# If the script is run without a command line argument, it regenerates a new
+# inputfile according to the initial conditions listed below.
+# For detailed references see:
+# 1) M. Tamburini. PhD Thesis. University of Pisa (2011)
+# https://etd.adm.unipi.it/t/etd-11232011-111059/
+# 2) H. Spohn Europhysics Letters 50 287 (2000)
+# https://arxiv.org/abs/physics/9911027
+# 3) H. Spohn, Dynamics of charged particles and their radiation field
+# (Cambridge University Press, Cambridge, 2004)
+
+import numpy as np
+import sys
+import yt
+
+#Input filename
+inputname = "inputs"
+#________________________________________
+
+#Physical constants
+c = 299792458.
+m_e = 9.1093837015e-31
+q_0 = 1.602176634e-19
+classical_electron_radius = 2.81794e-15
+reference_length = 1.0e-6
+very_small_dot_product = 1.0e-4
+very_small_weight = 1.0e-8
+#________________________________________
+
+#Sim box data
+sim_size = 0.8e-6
+resolution = 64
+steps = 64
+init_pos = np.array([0.,0.,0.])
+#________________________________________
+
+#Momentum vals
+p_aux_0 = np.array([2.,3.,6.])
+p_aux_1 = np.array([1,0,0])
+p_aux_2 = np.array([0,1,0])
+Q, _ = np.linalg.qr(np.column_stack( [p_aux_0, p_aux_1, p_aux_2] )) #Gram-Schmidt
+p_0 = -Q[:,0]
+p_1 = -Q[:,1]
+p_2 = -Q[:,2]
+
+p_vals = [50,200,1000]
+#________________________________________
+
+#Field val
+B_val_norm = 300
+B_val = B_val_norm*m_e * 2.0 * np.pi * c /q_0/reference_length
+B = p_0 * B_val
+#________________________________________
+
+#Tolerance
+tol = 0.05
+#________________________________________
+
+#tau_c
+omega_c = q_0*B_val/m_e
+t0 = (2./3.)*classical_electron_radius/c
+tau_c = 1.0/omega_c/omega_c/t0
+#________________________________________
+
+
+#Simulation case struct
+class sim_case:
+ def __init__(self, _name, _init_mom, _type):
+ self.name = _name
+ self.init_mom = _init_mom
+ self.type = _type
+#________________________________________
+
+#All cases
+cases = [
+ sim_case("ele_para0", p_0 * p_vals[2], "-q_e"),
+ sim_case("ele_perp0", p_1 * p_vals[0], "-q_e"),
+ sim_case("ele_perp1", p_2 * p_vals[1], "-q_e"),
+ sim_case("ele_perp2", p_1 * p_vals[2], "-q_e"),
+ sim_case("pos_perp2", p_1 * p_vals[2], "q_e"),
+]
+#________________________________________
+
+#Auxiliary functions
+def gamma(p) :
+ return np.sqrt(1.0 + np.dot(p,p))
+
+def exp_res(cc, time):
+ if np.all(np.linalg.norm(np.cross(cc.init_mom, B)) < very_small_dot_product):
+ return gamma(cc.init_mom)
+ else :
+ tt = time/tau_c
+ g0 = gamma(cc.init_mom)
+ C = -0.5 * np.log((g0+1)/(g0-1))
+ return 1.0/np.tanh(tt - C)
+#________________________________________
+
+
+def check():
+ filename = sys.argv[1]
+ data_set_end = yt.load(filename)
+
+ sim_time = data_set_end.current_time.to_value()
+
+ #simulation results
+ all_data = data_set_end.all_data()
+ spec_names = [cc.name for cc in cases]
+
+ #All momenta
+ res_mom = np.array([np.array([
+ all_data[sp, 'particle_momentum_x'].v[0],
+ all_data[sp, 'particle_momentum_y'].v[0],
+ all_data[sp, 'particle_momentum_z'].v[0]])
+ for sp in spec_names])
+
+ for cc in zip(cases, res_mom):
+ init_gamma = gamma(cc[0].init_mom)
+ end_gamma = gamma(cc[1]/m_e/c)
+ exp_gamma = exp_res(cc[0], sim_time)
+ assert(np.abs(end_gamma-exp_gamma)/exp_gamma < tol)
+
+def generate():
+
+ with open(inputname,'w') as f:
+ f.write("#Automatically generated inputfile\n")
+ f.write("#Run check.py without arguments to regenerate\n")
+ f.write("#\n\n")
+ f.write("max_step = {}\n".format(steps))
+ f.write("amr.n_cell = {} {} {}\n".format(resolution, resolution, resolution))
+ f.write("amr.max_level = 0\n")
+ f.write("amr.blocking_factor = 32\n")
+ f.write("amr.max_grid_size = 64\n")
+ f.write("geometry.coord_sys = 0\n")
+ f.write("geometry.is_periodic = 1 1 1\n")
+ f.write("geometry.prob_lo = {} {} {}\n".format(-sim_size, -sim_size, -sim_size))
+ f.write("geometry.prob_hi = {} {} {}\n".format(sim_size, sim_size, sim_size))
+ f.write("warpx.do_pml = 0\n")
+ f.write("algo.charge_deposition = standard\n")
+ f.write("algo.field_gathering = standard\n")
+ f.write("warpx.cfl = 1.0\n")
+ f.write("warpx.serialize_ics = 1\n")
+
+ f.write("\nparticles.nspecies = {}\n".format(len(cases)))
+ f.write("particles.species_names = ")
+ for cc in cases:
+ f.write(" {}".format(cc.name))
+ f.write("\n")
+
+ f.write("\namr.plot_int = {}\n\n".format(steps))
+
+ f.write("warpx.fields_to_plot = rho\n\n")
+
+ for cc in cases:
+ f.write("{}.charge = {}\n".format(cc.name, cc.type))
+ f.write("{}.mass = m_e\n".format(cc.name))
+ f.write('{}.injection_style = "SingleParticle"\n'.format(cc.name))
+ f.write("{}.single_particle_pos = {} {} {}\n".
+ format(cc.name, init_pos[0], init_pos[1], init_pos[2]))
+ f.write("{}.single_particle_vel = {} {} {}\n".
+ format(cc.name, cc.init_mom[0], cc.init_mom[1], cc.init_mom[2]))
+ f.write("{}.single_particle_weight = {}\n".format(cc.name, very_small_weight))
+ f.write("{}.do_classical_radiation_reaction = 1\n".format(cc.name))
+ f.write("\n")
+
+ f.write("warpx.B_external = {} {} {}\n".format(*B))
+
+
+def main():
+ if (len(sys.argv) < 2):
+ generate()
+ else:
+ check()
+
+if __name__ == "__main__":
+ main()
diff --git a/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs b/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs
new file mode 100644
index 000000000..51d89d740
--- /dev/null
+++ b/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs
@@ -0,0 +1,67 @@
+#Automatically generated inputfile
+#Run check.py without arguments to regenerate
+#
+
+max_step = 64
+amr.n_cell = 64 64 64
+amr.max_level = 0
+amr.blocking_factor = 32
+amr.max_grid_size = 64
+geometry.coord_sys = 0
+geometry.is_periodic = 1 1 1
+geometry.prob_lo = -8e-07 -8e-07 -8e-07
+geometry.prob_hi = 8e-07 8e-07 8e-07
+warpx.do_pml = 0
+algo.charge_deposition = standard
+algo.field_gathering = standard
+warpx.cfl = 1.0
+warpx.serialize_ics = 1
+
+particles.nspecies = 5
+particles.species_names = ele_para0 ele_perp0 ele_perp1 ele_perp2 pos_perp2
+
+amr.plot_int = 64
+
+warpx.fields_to_plot = rho
+
+ele_para0.charge = -q_e
+ele_para0.mass = m_e
+ele_para0.injection_style = "SingleParticle"
+ele_para0.single_particle_pos = 0.0 0.0 0.0
+ele_para0.single_particle_vel = 285.7142857142856 428.5714285714285 857.142857142857
+ele_para0.single_particle_weight = 1e-08
+ele_para0.do_classical_radiation_reaction = 1
+
+ele_perp0.charge = -q_e
+ele_perp0.mass = m_e
+ele_perp0.injection_style = "SingleParticle"
+ele_perp0.single_particle_pos = 0.0 0.0 0.0
+ele_perp0.single_particle_vel = -47.91574237499548 6.388765649999412 12.777531299998806
+ele_perp0.single_particle_weight = 1e-08
+ele_perp0.do_classical_radiation_reaction = 1
+
+ele_perp1.charge = -q_e
+ele_perp1.mass = m_e
+ele_perp1.injection_style = "SingleParticle"
+ele_perp1.single_particle_pos = 0.0 0.0 0.0
+ele_perp1.single_particle_vel = 1.4274296030894868e-14 178.88543819998318 -89.44271909999159
+ele_perp1.single_particle_weight = 1e-08
+ele_perp1.do_classical_radiation_reaction = 1
+
+ele_perp2.charge = -q_e
+ele_perp2.mass = m_e
+ele_perp2.injection_style = "SingleParticle"
+ele_perp2.single_particle_pos = 0.0 0.0 0.0
+ele_perp2.single_particle_vel = -958.3148474999095 127.77531299998823 255.55062599997612
+ele_perp2.single_particle_weight = 1e-08
+ele_perp2.do_classical_radiation_reaction = 1
+
+pos_perp2.charge = q_e
+pos_perp2.mass = m_e
+pos_perp2.injection_style = "SingleParticle"
+pos_perp2.single_particle_pos = 0.0 0.0 0.0
+pos_perp2.single_particle_vel = -958.3148474999095 127.77531299998823 255.55062599997612
+pos_perp2.single_particle_weight = 1e-08
+pos_perp2.do_classical_radiation_reaction = 1
+
+warpx.B_external = 917978.2333474257 1376967.350021139 2753934.700042278
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
index 90210a94e..a3521cf6a 100644
--- a/Regression/WarpX-tests.ini
+++ b/Regression/WarpX-tests.ini
@@ -650,6 +650,21 @@ doVis = 0
compareParticles = 0
analysisRoutine = Examples/Tests/photon_pusher/check.py
+
+[radiation_reaction]
+buildDir = .
+inputFile = Examples/Tests/radiation_reaction/test_const_B_analytical/inputs
+dim = 3
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 2
+compileTest = 0
+doVis = 0
+compareParticles = 0
+analysisRoutine = Examples/Tests/radiation_reaction/test_const_B_analytical/check.py
+
[qed_breit_wheeler_tau_init]
buildDir = .
inputFile = Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init
diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H
index 407cf26f3..9132cf4a9 100644
--- a/Source/Particles/PhotonParticleContainer.H
+++ b/Source/Particles/PhotonParticleContainer.H
@@ -73,6 +73,12 @@ public:
int lev,
int depos_lev,
amrex::Real dt) override {};
+ //Photons are not leptons
+ virtual bool AmIALepton () override
+ {
+ return false;
+ };
+
};
#endif // #ifndef WARPX_PhotonParticleContainer_H_
diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H
index 08493edcd..ed3a698e2 100644
--- a/Source/Particles/PhysicalParticleContainer.H
+++ b/Source/Particles/PhysicalParticleContainer.H
@@ -210,6 +210,14 @@ protected:
// Inject particles during the whole simulation
void ContinuousInjection (const amrex::RealBox& injection_box) override;
+ //This function return true if the PhysicalParticleContainer contains electrons
+ //or positrons, false otherwise
+ virtual bool AmIALepton ();
+
+ //When true PhysicalParticleContainer tries to use a pusher including
+ //radiation reaction
+ bool do_classical_radiation_reaction = false;
+
#ifdef WARPX_QED
// A flag to enable quantum_synchrotron process for leptons
bool do_qed_quantum_sync = false;
diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 66331db26..d12a4dbff 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -17,6 +17,7 @@
#include <UpdatePosition.H>
#include <UpdateMomentumBoris.H>
#include <UpdateMomentumVay.H>
+#include <UpdateMomentumBorisWithRadiationReaction.H>
#include <UpdateMomentumHigueraCary.H>
using namespace amrex;
@@ -45,6 +46,23 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp
pp.query("do_boosted_frame_diags", do_boosted_frame_diags);
pp.query("do_field_ionization", do_field_ionization);
+
+ //check if Radiation Reaction is enabled and do consistency checks
+ pp.query("do_classical_radiation_reaction", do_classical_radiation_reaction);
+ //if the species is not a lepton, do_classical_radiation_reaction
+ //should be false
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+ !(do_classical_radiation_reaction && !AmIALepton()),
+ "Can't enable classical radiation reaction for non lepton species. " );
+
+ //Only Boris pusher is compatible with radiation reaction
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
+ !(do_classical_radiation_reaction &&
+ WarpX::particle_pusher_algo != ParticlePusherAlgo::Boris),
+ "Radiation reaction can be enabled only if Boris pusher is used");
+ //_____________________________
+
+
#ifdef AMREX_USE_GPU
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(
do_field_ionization == 0,
@@ -1580,7 +1598,23 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti,
// Loop over the particles and update their momentum
const Real q = this->charge;
const Real m = this-> mass;
- if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
+
+
+ //Assumes that all consistency checks have been done at initialization
+ if(do_classical_radiation_reaction){
+ amrex::ParallelFor(
+ pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumBorisWithRadiationReaction( ux[i], uy[i], uz[i],
+ Ex[i], Ey[i], Ez[i], Bx[i],
+ By[i], Bz[i], qp, m, dt);
+ UpdatePosition( x[i], y[i], z[i],
+ ux[i], uy[i], uz[i], dt );
+ }
+ );
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
amrex::ParallelFor(
pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
@@ -1701,18 +1735,49 @@ PhysicalParticleContainer::PushP (int lev, Real dt,
// Loop over the particles and update their momentum
const Real q = this->charge;
const Real m = this-> mass;
- if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
+
+ int* AMREX_RESTRICT ion_lev = nullptr;
+ if (do_field_ionization){
+ ion_lev = pti.GetiAttribs(particle_icomps["ionization_level"]).dataPtr();
+ }
+
+
+ //Assumes that all consistency checks have been done at initialization
+ if(do_classical_radiation_reaction){
+ amrex::ParallelFor(
+ pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumBorisWithRadiationReaction(
+ ux[i], uy[i], uz[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ qp, m, dt);
+ }
+ );
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
amrex::ParallelFor( pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
- UpdateMomentumBoris( ux[i], uy[i], uz[i],
- Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumBoris(
+ ux[i], uy[i], uz[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ qp, m, dt);
}
);
} else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) {
amrex::ParallelFor( pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
- UpdateMomentumVay( ux[i], uy[i], uz[i],
- Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ Real qp = q;
+ if (ion_lev){ qp *= ion_lev[i]; }
+ UpdateMomentumVay(
+ ux[i], uy[i], uz[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ qp, m, dt);
}
);
} else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
@@ -2163,6 +2228,13 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const
);
}
+//This function return true if the PhysicalParticleContainer contains electrons
+//or positrons, false otherwise
+bool
+PhysicalParticleContainer::AmIALepton(){
+ return (this-> mass == PhysConst::m_e);
+}
+
#ifdef WARPX_QED
bool PhysicalParticleContainer::has_quantum_sync()
diff --git a/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
new file mode 100644
index 000000000..0abedd258
--- /dev/null
+++ b/Source/Particles/Pusher/UpdateMomentumBorisWithRadiationReaction.H
@@ -0,0 +1,84 @@
+
+#ifndef WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_
+#define WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_
+
+#include <AMReX_REAL.H>
+#include <UpdateMomentumBoris.H>
+
+/**
+ * Push the particle's positions over one timestep,
+ * given the value of its momenta `ux`, `uy`, `uz`.
+ * Includes Radiation Reaction according to
+ * https://doi.org/10.1088/1367-2630/12/12/123005
+ */
+AMREX_GPU_HOST_DEVICE AMREX_INLINE
+void UpdateMomentumBorisWithRadiationReaction(
+ amrex::ParticleReal& ux, amrex::ParticleReal& uy, amrex::ParticleReal& uz,
+ const amrex::ParticleReal Ex, const amrex::ParticleReal Ey, const amrex::ParticleReal Ez,
+ const amrex::ParticleReal Bx, const amrex::ParticleReal By, const amrex::ParticleReal Bz,
+ const amrex::Real q, const amrex::Real m, const amrex::Real dt )
+{
+ //RR algorithm needs to store old value of the normalized momentum
+ const amrex::Real ux_old = ux;
+ const amrex::Real uy_old = uy;
+ const amrex::Real uz_old = uz;
+
+ //Useful constants
+ const amrex::Real inv_dt = 1.0/dt;
+ constexpr amrex::Real inv_c2 = 1./(PhysConst::c*PhysConst::c);
+
+ //Call to regular Boris pusher
+ UpdateMomentumBoris(
+ ux, uy, uz,
+ Ex, Ey, Ez,
+ Bx, By, Bz,
+ q, m, dt );
+
+ //Estimation of the normalized momentum at intermediate (integer) time
+ const amrex::Real ux_n = (ux+ux_old)*0.5;
+ const amrex::Real uy_n = (uy+uy_old)*0.5;
+ const amrex::Real uz_n = (uz+uz_old)*0.5;
+
+ // Compute Lorentz factor (and inverse) at intermediate (integer) time
+ const amrex::Real gamma_n = std::sqrt( 1. +
+ (ux_n*ux_n + uy_n*uy_n + uz_n*uz_n)*inv_c2);
+ const amrex::Real inv_gamma_n = 1.0/gamma_n;
+
+ //Estimation of the velocity at intermediate (integer) time
+ const amrex::Real vx_n = ux_n*inv_gamma_n;
+ const amrex::Real vy_n = uy_n*inv_gamma_n;
+ const amrex::Real vz_n = uz_n*inv_gamma_n;
+ const amrex::Real bx_n = vx_n/PhysConst::c;
+ const amrex::Real by_n = vy_n/PhysConst::c;
+ const amrex::Real bz_n = vz_n/PhysConst::c;
+
+ //Lorentz force over charge
+ const amrex::Real flx_q = (Ex + vy_n*Bz - vz_n*By);
+ const amrex::Real fly_q = (Ey + vz_n*Bx - vx_n*Bz);
+ const amrex::Real flz_q = (Ez + vx_n*By - vy_n*Bx);
+ const amrex::Real fl_q2 = flx_q*flx_q + fly_q*fly_q + flz_q*flz_q;
+
+ //Calculation of auxiliary quantities
+ const amrex::Real bdotE = (bx_n*Ex + by_n*Ey + bz_n*Ez);
+ const amrex::Real bdotE2 = bdotE*bdotE;
+ const amrex::Real coeff = gamma_n*gamma_n*(fl_q2-bdotE2);
+
+ //Radiation reaction constant
+ const amrex::Real RRcoeff = 2.0*PhysConst::r_e*q*q/
+ (3.0*m*m*PhysConst::c*PhysConst::c);
+
+ //Compute the components of the RR force
+ const amrex::Real frx =
+ RRcoeff*(PhysConst::c*(fly_q*Bz - flz_q*By) + bdotE*Ex - coeff*bx_n);
+ const amrex::Real fry =
+ RRcoeff*(PhysConst::c*(flz_q*Bx - flx_q*Bz) + bdotE*Ey - coeff*by_n);
+ const amrex::Real frz =
+ RRcoeff*(PhysConst::c*(flx_q*By - fly_q*Bx) + bdotE*Ez - coeff*bz_n);
+
+ //Update momentum using the RR force
+ ux += frx*dt;
+ uy += fry*dt;
+ uz += frz*dt;
+}
+
+#endif // WARPX_PARTICLES_PUSHER_UPDATEMOMENTUM_BORIS_WITHRR_H_
diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp
index 2ef833151..535ffec6f 100644
--- a/Source/Particles/RigidInjectedParticleContainer.cpp
+++ b/Source/Particles/RigidInjectedParticleContainer.cpp
@@ -13,6 +13,7 @@
#include <WarpXAlgorithmSelection.H>
#include <UpdateMomentumBoris.H>
#include <UpdateMomentumVay.H>
+#include <UpdateMomentumBorisWithRadiationReaction.H>
#include <UpdateMomentumHigueraCary.H>
using namespace amrex;
@@ -430,18 +431,37 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt,
// Loop over the particles and update their momentum
const Real q = this->charge;
const Real m = this->mass;
- if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
+
+ //Assumes that all consistency checks have been done at initialization
+ if(do_classical_radiation_reaction){
+ amrex::ParallelFor(
+ pti.numParticles(),
+ [=] AMREX_GPU_DEVICE (long i) {
+ UpdateMomentumBorisWithRadiationReaction(
+ uxpp[i], uypp[i], uzpp[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ q, m, dt);
+ }
+ );
+ } else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){
amrex::ParallelFor( pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
- UpdateMomentumBoris( uxpp[i], uypp[i], uzpp[i],
- Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ UpdateMomentumBoris(
+ uxpp[i], uypp[i], uzpp[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ q, m, dt);
}
);
} else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Vay) {
amrex::ParallelFor( pti.numParticles(),
[=] AMREX_GPU_DEVICE (long i) {
- UpdateMomentumVay( uxpp[i], uypp[i], uzpp[i],
- Expp[i], Eypp[i], Ezpp[i], Bxpp[i], Bypp[i], Bzpp[i], q, m, dt);
+ UpdateMomentumVay(
+ uxpp[i], uypp[i], uzpp[i],
+ Expp[i], Eypp[i], Ezpp[i],
+ Bxpp[i], Bypp[i], Bzpp[i],
+ q, m, dt);
}
);
} else if (WarpX::particle_pusher_algo == ParticlePusherAlgo::HigueraCary) {
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-19 06:43:54 +0000