Hi-PACE · SeverinDiederichs · Jun 14, 2023 · Jun 8, 2023 · Jun 9, 2023 · Jun 9, 2023
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -257,6 +257,12 @@ if(BUILD_TESTING)
                 WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
         )
 
+        add_test(NAME radiation_reactions.1Rank
+                COMMAND bash ${HiPACE_SOURCE_DIR}/tests/radiation_reactions.1Rank.sh
+                        $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}
+                WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
+        )
+
         add_test(NAME grid_current.1Rank
                 COMMAND bash ${HiPACE_SOURCE_DIR}/tests/grid_current.1Rank.sh
                         $<TARGET_FILE:HiPACE> ${HiPACE_SOURCE_DIR}

diff --git a/docs/source/run/parameters.rst b/docs/source/run/parameters.rst
@@ -518,6 +518,15 @@ which are valid only for certain beam types, are introduced further below under
     For the last component, z actually represents the zeta coordinate zeta = z - c*t.
     For instance, ``hipace.external_B_slope = -1. 1. 0.`` creates an axisymmetric focusing lens of strength 1 T/m.
 
+* ``<beam name>.do_z_push`` (`bool`) optional (default `1`)
+    Whether the beam particles are pushed along the z-axis. The momentum is still fully updated.
+    Note: using ``do_z_push = 0`` results in unphysical behavior.
+
+* ``<beam name> or beams..do_radiation_reaction`` (`bool`) optional (default `0`)
+    Whether the beam particles undergo energy loss due to classical radiation reactions.
+    The implemented radiation reaction model is based on this publication: `M. Tamburini et al., NJP 12, 123005 <https://doi.org/10.1088/1367-2630/12/12/123005>`__
+    Currently only implemented for SI units, normalized units will be added soon.
+
 Option: ``fixed_weight``
 ^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -546,10 +555,6 @@ Option: ``fixed_weight``
     ``beam_name.num_particles``, therefore this option requires that the beam particle number must be
     divisible by 4.
 
-* ``<beam name>.do_z_push`` (`bool`) optional (default `1`)
-    Whether the beam particles are pushed along the z-axis. The momentum is still fully updated.
-    Note: using ``do_z_push = 0`` results in unphysical behavior.
-
 * ``<beam name>.z_foc`` (`float`) optional (default `0.`)
     Distance at which the beam will be focused, calculated from the position at which the beam is initialized.
     The beam is assumed to propagate ballistically in-between.

diff --git a/examples/beam_in_vacuum/analysis_RR.py b/examples/beam_in_vacuum/analysis_RR.py
@@ -0,0 +1,62 @@
+#! /usr/bin/env python3
+
+# Copyright 2020-2023
+#
+# This file is part of HiPACE++. It tests the radiation reaction of a beam in an external field vs
+# the analytic theory in P. Michel et al., PRE 74, 026501 https://doi.org/10.1103/PhysRevE.74.026501
+#
+# Authors: Severin Diederichs
+# License: BSD-3-Clause-LBNL
+
+
+import numpy as np
+import sys
+sys.path.append("../../tools/")
+import read_insitu_diagnostics as diag
+from scipy import constants as scc
+
+# load HiPACE++ data with insitu diags
+all_data_with_RR = diag.read_file('diags/insitu/reduced_beam.*.txt')
+
+ne = 5e24
+wp = np.sqrt(ne*scc.e**2 / (scc.m_e * scc.epsilon_0 ))
+kp = wp/scc.c
+
+mean_gamma0 = diag.gamma_mean(all_data_with_RR["average"])[0] # should be 2000
+std_gamma0 = diag.gamma_spread(all_data_with_RR["average"])[0] \
+            /diag.gamma_mean(all_data_with_RR["average"])[0] # should be 0.01
+epsilonx0 = diag.emittance_x(all_data_with_RR["average"])[0] # sigma_x0**2 *np.sqrt(mean_gamma0/2) *kp
+# should be 313e-6
+
+# final simulation values
+mean_gamma_sim = diag.gamma_mean(all_data_with_RR["average"])[-1]
+std_gamma_sim = diag.gamma_spread(all_data_with_RR["average"])[-1] \
+            /diag.gamma_mean(all_data_with_RR["average"])[-1]
+epsilonx_sim = diag.emittance_x(all_data_with_RR["average"])[-1]
+
+# calculate theoretical values
+sigma_x0 = np.sqrt(epsilonx0 / (kp* np.sqrt(mean_gamma0/2))) # should be 4.86e-6
+ux0 = epsilonx0/sigma_x0
+r_e = 1/(4*np.pi*scc.epsilon_0) * scc.e**2/(scc.m_e*scc.c**2)
+taur = 6.24e-24 # 2*r_e /(3*scc.c)
+K = kp/np.sqrt(2)
+w_beta = K*scc.c/np.sqrt(mean_gamma0)
+xmsq = sigma_x0**2 + scc.c**2*ux0**2/(w_beta**2 * mean_gamma0**2)
+nugamma = taur * scc.c**2 * K**4 * mean_gamma0 * xmsq/2 # equation 32 from the paper
+nugammastd = taur * scc.c**2 * K**4 * mean_gamma0 * sigma_x0**2
+
+t = all_data_with_RR["time"][-1]
+gamma_theo = mean_gamma0/(1+nugamma*t) # equation 31 from the paper
+std_gamma_theo = np.sqrt(std_gamma0**2 + nugammastd**2 * t**2) # equation 35 from the paper
+emittance_theo = epsilonx0/(1+3*nugammastd*t/2) # equation 39 from the paper
+
+error_analytic_gamma = np.abs(mean_gamma_sim-gamma_theo)/gamma_theo
+error_analytic_std_gamma = np.abs(std_gamma_sim-std_gamma_theo)/std_gamma_theo
+error_analytic_emittance = np.abs(epsilonx_sim-emittance_theo)/emittance_theo
+
+assert(error_analytic_gamma < 1e-3)
+assert(error_analytic_std_gamma < 3e-2)
+assert(error_analytic_emittance < 1e-3)
+print("Error on gamma ", error_analytic_gamma)
+print("Error on relative gamma spread ", error_analytic_std_gamma)
+print("Error on emittance ", error_analytic_emittance)
diff --git a/examples/beam_in_vacuum/inputs_RR b/examples/beam_in_vacuum/inputs_RR
@@ -0,0 +1,43 @@
+amr.n_cell = 16 16 4
+
+amr.max_level = 0
+
+my_constants.ne = 5e24
+my_constants.wp = sqrt( ne * q_e^2 / (epsilon0 * m_e))
+my_constants.E0 = wp * m_e * clight / q_e
+my_constants.kp = wp / clight
+my_constants.kp_inv = clight / wp
+
+my_constants.K = kp/sqrt(2.)
+my_constants.gamma0 = 2000
+my_constants.emittance_x = 313e-6
+my_constants.sigma_x = sqrt(emittance_x*kp_inv / sqrt(gamma0/2.) )
+my_constants.sigma_ux = emittance_x / sigma_x
+my_constants.uz = sqrt(gamma0^2 - 1 - sigma_ux^2)
+my_constants.w_beta = K*clight/sqrt(gamma0)
+
+hipace.external_E_slope = 1/2*E0/kp_inv 1/2*E0/kp_inv 0.
+
+hipace.dt = 30 /w_beta
+hipace.verbose = 1
+max_step = 5
+diagnostic.output_period = 5
+diagnostic.diag_type = xz
+
+geometry.is_periodic = 1     1     1      # Is periodic?
+geometry.prob_lo     = -30.e-6   -30.e-6   -10.e-6    # physical domain
+geometry.prob_hi     =  30.e-6    30.e-6    10.e-6
+
+beams.names = beam
+beams.insitu_period = 1
+beam.injection_type = fixed_weight
+beam.profile = gaussian
+beam.position_mean = 0 0 0
+beam.position_std = sigma_x 1e-12 1e-6
+beam.density = ne/1e10
+beam.u_mean = 0. 0. uz
+beam.u_std = sigma_ux 0 uz*0.01
+beam.num_particles = 100000
+beam.n_subcycles = 50
+beam.do_radiation_reaction = 1
+beam.do_z_push = 0 # to avoid dephasing
diff --git a/src/particles/beam/BeamParticleContainer.H b/src/particles/beam/BeamParticleContainer.H
@@ -159,6 +159,7 @@ public:
     amrex::Real m_mass; /**< mass of each particle of this species */
     bool m_do_z_push {true}; /**< Pushing beam particles in z direction */
     int m_n_subcycles {10}; /**< Number of sub-cycles in the beam pusher */
+    bool m_do_radiation_reaction {false}; /**< whether to calculate radiation losses */
     /** Number of particles on upstream rank (required for IO) */
     bool m_do_salame = false; /**< Whether this beam uses salame */
     int m_num_particles_on_upstream_ranks {0};

diff --git a/src/particles/beam/BeamParticleContainer.cpp b/src/particles/beam/BeamParticleContainer.cpp
@@ -57,6 +57,10 @@ BeamParticleContainer::ReadParameters ()
     getWithParser(pp, "injection_type", m_injection_type);
     queryWithParser(pp, "duz_per_uz0_dzeta", m_duz_per_uz0_dzeta);
     queryWithParser(pp, "do_z_push", m_do_z_push);
+    queryWithParserAlt(pp, "do_radiation_reaction", m_do_radiation_reaction, pp_alt);
+    if (m_do_radiation_reaction) AMREX_ALWAYS_ASSERT_WITH_MESSAGE( Hipace::m_normalized_units == 0,
+        "Radiation reactions are only implemented in SI units so far. Please use  "
+        "`hipace.normalized_units = 0`");
     queryWithParserAlt(pp, "insitu_period", m_insitu_period, pp_alt);
     queryWithParserAlt(pp, "insitu_file_prefix", m_insitu_file_prefix, pp_alt);
     queryWithParser(pp, "n_subcycles", m_n_subcycles);

diff --git a/src/particles/pusher/BeamParticleAdvance.cpp b/src/particles/pusher/BeamParticleAdvance.cpp
@@ -28,6 +28,7 @@ AdvanceBeamParticlesSlice (
 
     const bool do_z_push = beam.m_do_z_push;
     const int n_subcycles = beam.m_n_subcycles;
+    const bool radiation_reaction = beam.m_do_radiation_reaction;
     const amrex::Real dt = Hipace::m_dt / n_subcycles;
 
     const int psi_comp = Comps[WhichSlice::This]["Psi"];
@@ -98,7 +99,8 @@ AdvanceBeamParticlesSlice (
     int const num_particles = cell_stop-cell_start;
 
     const amrex::Real clight = phys_const.c;
-    const amrex::Real clightsq = 1.0_rt/(phys_const.c*phys_const.c);
+    const amrex::Real inv_clight = 1.0_rt/phys_const.c;
+    const amrex::Real inv_c2 = 1.0_rt/(phys_const.c*phys_const.c);
     const amrex::Real charge_mass_ratio = beam.m_charge / beam.m_mass;
     const amrex::RealVect external_E_uniform = extEu;
     const amrex::RealVect external_B_uniform = extBu;
@@ -124,9 +126,7 @@ AdvanceBeamParticlesSlice (
             for (int i = 0; i < n_subcycles; i++) {
 
                 const amrex::ParticleReal gammap_inv = 1._rt / std::sqrt( 1._rt
-                    + ux*ux*clightsq
-                    + uy*uy*clightsq
-                    + uz*uz*clightsq );
+                    + (ux*ux + uy*uy + uz*uz)*inv_c2 );
 
                 // first we do half a step in x,y
                 // This is not required in z, which is pushed in one step later
@@ -174,9 +174,9 @@ AdvanceBeamParticlesSlice (
                     external_E_uniform, external_B_uniform, external_E_slope, external_B_slope);
 
                 // use intermediate fields to calculate next (n+1) transverse momenta
-                const amrex::ParticleReal ux_next = ux + dt * charge_mass_ratio
+                amrex::ParticleReal ux_next = ux + dt * charge_mass_ratio
                     * ( ExmByp + ( clight - uz * gammap_inv ) * Byp + uy*gammap_inv*Bzp);
-                const amrex::ParticleReal uy_next = uy + dt * charge_mass_ratio
+                amrex::ParticleReal uy_next = uy + dt * charge_mass_ratio
                     * ( EypBxp + ( uz * gammap_inv - clight ) * Bxp - ux*gammap_inv*Bzp);
 
                 // Now computing new longitudinal momentum
@@ -186,19 +186,66 @@ AdvanceBeamParticlesSlice (
                     + dt * 0.5_rt * charge_mass_ratio * Ezp;
 
                 const amrex::ParticleReal gamma_intermediate_inv = 1._rt / std::sqrt( 1._rt
-                    + ux_intermediate*ux_intermediate*clightsq
-                    + uy_intermediate*uy_intermediate*clightsq
-                    + uz_intermediate*uz_intermediate*clightsq );
+                    + (ux_intermediate*ux_intermediate
+                    +  uy_intermediate*uy_intermediate
+                    +  uz_intermediate*uz_intermediate)*inv_c2 );
 
-                const amrex::ParticleReal uz_next = uz + dt * charge_mass_ratio
+                amrex::ParticleReal uz_next = uz + dt * charge_mass_ratio
                     * ( Ezp + ( ux_intermediate * Byp - uy_intermediate * Bxp )
                     * gamma_intermediate_inv );
 
+                if (radiation_reaction) {
+
+                    const amrex::ParticleReal Exp = ExmByp + clight*Byp;
+                    const amrex::ParticleReal Eyp = EypBxp - clight*Bxp;
+
+                    const amrex::ParticleReal gamma_intermediate = std::sqrt( 1._rt
+                        + (ux_intermediate*ux_intermediate
+                        +  uy_intermediate*uy_intermediate
+                        +  uz_intermediate*uz_intermediate)*inv_c2 );
+                    // Estimation of the velocity at intermediate time
+                    const amrex::ParticleReal vx_n = ux_intermediate*gamma_intermediate_inv;
+                    const amrex::ParticleReal vy_n = uy_intermediate*gamma_intermediate_inv;
+                    const amrex::ParticleReal vz_n = uz_intermediate*gamma_intermediate_inv;
+                    const amrex::ParticleReal bx_n = vx_n*inv_clight;
+                    const amrex::ParticleReal by_n = vy_n*inv_clight;
+                    const amrex::ParticleReal bz_n = vz_n*inv_clight;
+
+                    // Lorentz force over charge
+                    const amrex::ParticleReal flx_q = (Exp + vy_n*Bzp - vz_n*Byp);
+                    const amrex::ParticleReal fly_q = (Eyp + vz_n*Bxp - vx_n*Bzp);
+                    const amrex::ParticleReal flz_q = (Ezp + vx_n*Byp - vy_n*Bxp);
+                    const amrex::ParticleReal fl_q2 = flx_q*flx_q + fly_q*fly_q + flz_q*flz_q;
+
+                    // Calculation of auxiliary quantities
+                    const amrex::ParticleReal bdotE = (bx_n*Exp + by_n*Eyp + bz_n*Ezp);
+                    const amrex::ParticleReal bdotE2 = bdotE*bdotE;
+                    const amrex::ParticleReal coeff = gamma_intermediate*gamma_intermediate*(fl_q2-bdotE2);
+
+                    // Radiation reaction constant
+                    const amrex::ParticleReal q_over_mc = charge_mass_ratio*inv_clight;
+                    const amrex::ParticleReal RRcoeff = (2.0_rt/3.0_rt)*PhysConstSI::r_e*q_over_mc*q_over_mc;
+
+                    //Compute the components of the RR force
+                    const amrex::ParticleReal frx =
+                        RRcoeff*(PhysConstSI::c*(fly_q*Bzp - flz_q*Byp) + bdotE*Exp - coeff*bx_n);
+                    const amrex::ParticleReal fry =
+                        RRcoeff*(PhysConstSI::c*(flz_q*Bxp - flx_q*Bzp) + bdotE*Eyp - coeff*by_n);
+                    const amrex::ParticleReal frz =
+                        RRcoeff*(PhysConstSI::c*(flx_q*Byp - fly_q*Bxp) + bdotE*Ezp - coeff*bz_n);
+
+                    //Update momentum using the RR force
+                    ux_next += frx*dt;
+                    uy_next += fry*dt;
+                    uz_next += frz*dt;
+                }
+
+
                 /* computing next gamma value */
                 const amrex::ParticleReal gamma_next_inv = 1._rt / std::sqrt( 1._rt
-                    + ux_next*ux_next*clightsq
-                    + uy_next*uy_next*clightsq
-                    + uz_next*uz_next*clightsq );
+                    + (ux_next*ux_next
+                    +  uy_next*uy_next
+                    +  uz_next*uz_next)*inv_c2 );
 
                 /*
                  * computing positions and setting momenta for the next timestep

diff --git a/src/utils/Constants.H b/src/utils/Constants.H
@@ -21,6 +21,7 @@ namespace PhysConstSI
     static constexpr auto m_e = static_cast<amrex::Real>( 9.1093837015e-31 );
     static constexpr auto m_p = static_cast<amrex::Real>( 1.67262192369e-27 );
     static constexpr auto hbar  = static_cast<amrex::Real>( 1.054571817e-34 );
+    static constexpr auto r_e = static_cast<amrex::Real>( 2.817940326204929e-15 );
 }
 
 /** \brief Namespace containing math constants */

diff --git a/tests/adaptive_time_step.1Rank.sh b/tests/adaptive_time_step.1Rank.sh
@@ -22,11 +22,14 @@ HIPACE_SOURCE_DIR=$2
 HIPACE_EXAMPLE_DIR=${HIPACE_SOURCE_DIR}/examples/beam_in_vacuum
 HIPACE_TEST_DIR=${HIPACE_SOURCE_DIR}/tests
 
+FILE_NAME=`basename "$0"`
+TEST_NAME="${FILE_NAME%.*}"
+
 # Relative tolerance for checksum tests depends on the platform
 RTOL=1e-12 && [[ "$HIPACE_EXECUTABLE" == *"hipace"*".CUDA."* ]] && RTOL=1e-5
 
-rm -rf negative_gradient_data
-rm -rf positive_gradient_data
+rm -rf negative_gradient.txt
+rm -rf positive_gradient.txt
 
 # Run the simulation
 mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
@@ -36,7 +39,7 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
         geometry.prob_lo = -2. -2. -2. \
         geometry.prob_hi =  2.  2.  2. \
         hipace.dt = 0. \
-        diagnostic.output_period = 20 \
+        diagnostic.output_period = 0 \
         beam.density = 1 \
         beam.radius = 1. \
         beam.n_subcycles = 4 \
@@ -65,6 +68,7 @@ mpiexec -n 1 $HIPACE_EXECUTABLE $HIPACE_EXAMPLE_DIR/inputs_normalized \
         hipace.dt = adaptive\
         plasmas.adaptive_density=1 \
         hipace.nt_per_betatron = 89.7597901025655 \
+        hipace.file_prefix=$TEST_NAME \
         > positive_gradient.txt
 
 # Compare the result with theory
@@ -74,5 +78,5 @@ $HIPACE_EXAMPLE_DIR/analysis_adaptive_ts.py
 $HIPACE_TEST_DIR/checksum/checksumAPI.py \
     --evaluate \
     --rtol $RTOL \
-    --file_name diags/hdf5 \
-    --test-name adaptive_time_step.1Rank
+    --file_name $TEST_NAME \
+    --test-name $TEST_NAME
diff --git a/tests/checksum/benchmarks_json/radiation_reactions.1Rank.json b/tests/checksum/benchmarks_json/radiation_reactions.1Rank.json
@@ -0,0 +1,32 @@
+{
+  "lev=0": {
+    "Bx": 3.9728499837158e-19,
+    "By": 5.8910609879039e-14,
+    "Bz": 1.6891286774736e-24,
+    "ExmBy": 0.0,
+    "EypBx": 0.0,
+    "Ez": 3.1259472441681e-09,
+    "Psi": 0.0,
+    "Sx": 0.0024041543676677,
+    "Sy": 7.6433815475634e-12,
+    "chi": 0.0,
+    "jx": 3.1950444198537e-06,
+    "jx_beam": 3.1950444198537e-06,
+    "jy": 3.7374939092852e-13,
+    "jy_beam": 3.7374939092852e-13,
+    "jz_beam": 0.013037977938616,
+    "rhomjz": 0.0
+  },
+  "beam": {
+    "charge": 1.602176634e-14,
+    "id": 5000050000,
+    "mass": 9.1093837015e-26,
+    "x": 0.37632823449197,
+    "y": 4.6797140205514e-08,
+    "z": 0.079644759023598,
+    "ux": 5191150.7069722,
+    "uy": 0.7570356636811,
+    "uz": 196732680.56282,
+    "w": 3.8193025298947e-08
+  }
+}
diff --git a/tests/checksum/reset_all_benchmarks.sh b/tests/checksum/reset_all_benchmarks.sh
@@ -231,10 +231,22 @@ then
         || echo "ctest command failed, maybe just because checksums are different. Keep going"
     cd $checksum_dir
     ./checksumAPI.py --reset-benchmark \
-                     --file_name ${build_dir}/bin/diags/hdf5 \
+                     --file_name ${build_dir}/bin/adaptive_time_step.1Rank \
                      --test-name adaptive_time_step.1Rank
 fi
 
+# radiation_reactions.1Rank
+if [[ $all_tests = true ]] || [[ $one_test_name = "radiation_reactions.1Rank" ]]
+then
+    cd $build_dir
+    ctest --output-on-failure -R radiation_reactions.1Rank \
+        || echo "ctest command failed, maybe just because checksums are different. Keep going"
+    cd $checksum_dir
+    ./checksumAPI.py --reset-benchmark \
+                     --file_name ${build_dir}/bin/radiation_reactions.1Rank \
+                     --test-name radiation_reactions.1Rank
+fi
+
 # grid_current.1Rank
 if [[ $all_tests = true ]] || [[ $one_test_name = "grid_current.1Rank" ]]
 then