add classical radiation reactions for beams

CMakeLists.txt

@@ -257,6 +257,12 @@ if(BUILD_TESTING)
                 WORKING_DIRECTORY ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}
         )
 
+        add_test(NAME radiation_reaction.1Rank
+                COMMAND bash ${HiPACE_SOURCE_DIR}/tests/radiation_reaction.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}

docs/source/run/parameters.rst

@@ -218,6 +218,10 @@ General parameters
     The default value of this function corresponds to a flat Ez field at the position of the SALAME beam.
     Note: `zeta` is always less than or equal to `zeta_initial` and `Ez_initial` is typically below zero for electron beams.
 
+* ``hipace.background_density_SI`` (`float`) optional
+    Background plasma density in SI units. Certain physical modules (collisions, ionization, radiation reactions) depend on the actual background density.
+    Hence, in normalized units, they can only be included, if a background plasma density in SI units is provided using this input parameter.
+
 Field solver parameters
 -----------------------
 
@@ -412,10 +416,7 @@ Binary collisions for plasma species
 WARNING: this module is in development.
 
 HiPACE++ proposes an implementation of [Perez et al., Phys. Plasmas 19, 083104 (2012)], inherited from WarpX, between plasma species.
-
-* ``plasmas.background_density_SI`` (`float`) optional
-    Background plasma density in SI units. Only used for collisions in normalized units. Since the collision rate depends on the plasma density itself, it cannot be determined in normalized units without knowing the actual plasma background density.
-    Hence, it must be provided using this input parameter.
+As collisions depend on the physical density, in normalized units `hipace.background_density_SI` must be specified.
 
 * ``plasmas.collisions`` (list of `strings`) optional
     List of names of types binary Coulomb collisions.
@@ -518,6 +519,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>`__
+    In normalized units, `hipace.background_density_SI` must be specified.
+
 Option: ``fixed_weight``
 ^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -546,10 +556,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.

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
+
+print("Error on gamma ", error_analytic_gamma)
+print("Error on relative gamma spread ", error_analytic_std_gamma)
+print("Error on emittance ", error_analytic_emittance)
+assert(error_analytic_gamma < 1e-3)
+assert(error_analytic_std_gamma < 3e-2)
+assert(error_analytic_emittance < 1e-3)

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 = 1 / kp
+
+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*kp*E0 1/2*kp*E0 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

src/Hipace.H

@@ -335,6 +335,9 @@ public:
     inline static bool m_use_amrex_mlmg = false;
     /** Whether the simulation uses a laser pulse */
     inline static bool m_use_laser = false;
+    /** Background plasma density in SI, used to compute collisions, ionization,
+     * or radiation reaction in normalized units */
+    inline static amrex::Real m_background_density_SI = 0.;
     /** Adaptive time step instance */
     AdaptiveTimeStep m_adaptive_time_step;
     /** Laser instance (soon to be multi laser container) */

src/Hipace.cpp

@@ -162,6 +162,8 @@ Hipace::Hipace () :
     AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_do_tiling==0, "Tiling must be turned off to run on GPU.");
 #endif
 
+    queryWithParser(pph, "background_density_SI", m_background_density_SI);
+
 #ifdef AMREX_USE_MPI
     queryWithParser(pph, "skip_empty_comms", m_skip_empty_comms);
     int myproc = amrex::ParallelDescriptor::MyProc();

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};

src/particles/beam/BeamParticleContainer.cpp

@@ -57,6 +57,7 @@ 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);
     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);

src/particles/beam/MultiBeam.H

@@ -73,7 +73,7 @@ public:
                   The components represent d(Bx)/dy, d(By)/dx and d(Bz)/dz respectively.
                   Note the order of derivatives for the transverse components!
                   For the last component, z represents the zeta coordinate zeta = z - c*t
-    */
+     */
     void AdvanceBeamParticlesSlice (
         const Fields& fields, amrex::Vector<amrex::Geometry> const& gm, int const current_N_level,
         const int islice, const amrex::RealVect& extEu, const amrex::RealVect& extBu,

src/particles/plasma/MultiPlasma.H

@@ -168,8 +168,6 @@ private:
     std::vector<std::string> m_collision_names;
     /** Vector of binary collisions */
     amrex::Vector< CoulombCollision > m_all_collisions;
-    /** Background density in SI, used to compute collisions in normalized units */
-    amrex::Real m_background_density_SI = 0.;
 };
 
 #endif // MULTIPLASMA_H_

src/particles/plasma/MultiPlasma.cpp

@@ -11,6 +11,7 @@
 #include "particles/pusher/PlasmaParticleAdvance.H"
 #include "particles/sorting/TileSort.H"
 #include "utils/HipaceProfilerWrapper.H"
+#include "utils/DeprecatedInput.H"
 #include "Hipace.H"
 
 MultiPlasma::MultiPlasma ()
@@ -20,7 +21,9 @@ MultiPlasma::MultiPlasma ()
     queryWithParser(pp, "adaptive_density", m_adaptive_density);
     queryWithParser(pp, "sort_bin_size", m_sort_bin_size);
     queryWithParser(pp, "collisions", m_collision_names);
-    queryWithParser(pp, "background_density_SI", m_background_density_SI);
+
+    DeprecatedInput ("plasmas", "background_density_SI",
+                    "hipace.background_density_SI", "", true);
 
     if (m_names[0] == "no_plasma") return;
     m_nplasmas = m_names.size();
@@ -34,10 +37,10 @@ MultiPlasma::MultiPlasma ()
      for (int i = 0; i < m_ncollisions; ++i) {
          m_all_collisions.emplace_back(CoulombCollision(m_names, m_collision_names[i]));
      }
-     if (Hipace::GetInstance().m_normalized_units && m_ncollisions > 0) {
-         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_background_density_SI!=0,
+     if (Hipace::m_normalized_units && m_ncollisions > 0) {
+         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(Hipace::m_background_density_SI!=0,
              "For collisions with normalized units, a background plasma density must "
-             "be specified via 'plasmas.background_density_SI'");
+             "be specified via 'hipace.background_density_SI'");
      }
 }
 
@@ -61,7 +64,8 @@ MultiPlasma::InitData (amrex::Vector<amrex::BoxArray> slice_ba,
             }
             AMREX_ALWAYS_ASSERT_WITH_MESSAGE(plasma_product != nullptr,
                 "Must specify a valid product plasma for ionization using ionization_product");
-            plasma.InitIonizationModule(gm[0], plasma_product, m_background_density_SI); // geometry only for dz
+            plasma.InitIonizationModule(gm[0], plasma_product,
+            Hipace::m_background_density_SI); // geometry only for dz
         }
     }
     if (m_nplasmas > 0) m_all_bins.resize(m_nplasmas);
@@ -129,7 +133,7 @@ MultiPlasma::DoFieldIonization (
     const int lev, const amrex::Geometry& geom, const Fields& fields)
 {
     for (auto& plasma : m_all_plasmas) {
-        plasma.IonizationModule(lev, geom, fields, m_background_density_SI);
+        plasma.IonizationModule(lev, geom, fields, Hipace::m_background_density_SI);
     }
 }
 
@@ -177,7 +181,7 @@ MultiPlasma::doCoulombCollision (int lev, amrex::Box bx, amrex::Geometry geom)
 
         CoulombCollision::doCoulombCollision(
             lev, bx, geom, species1, species2, m_all_collisions[i].m_isSameSpecies,
-            m_all_collisions[i].m_CoulombLog, m_background_density_SI);
+            m_all_collisions[i].m_CoulombLog, Hipace::m_background_density_SI);
     }
 }
 

src/particles/plasma/PlasmaParticleContainerInit.cpp

@@ -326,7 +326,7 @@ InitIonizationModule (const amrex::Geometry& geom, PlasmaParticleContainer* prod
     if (normalized_units) {
         AMREX_ALWAYS_ASSERT_WITH_MESSAGE(background_density_SI!=0,
             "For ionization with normalized units, a background plasma density != 0 must "
-            "be specified via 'plasmas.background_density_SI'");
+            "be specified via 'hipace.background_density_SI'");
     }
 
     m_product_pc = product_pc;