From c1a553681da462cfe1fb797849e4bd3516e01866 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 3 Jan 2025 15:02:01 -0800 Subject: [PATCH 01/38] Initial implementation of Bremmstrahlung collisions --- .../BinaryCollision/BinaryCollision.H | 21 +- .../BinaryCollision/BinaryCollisionUtils.H | 1 + .../BinaryCollision/BinaryCollisionUtils.cpp | 3 + .../Bremsstrahlung/BremsstrahlungFunc.H | 477 ++++++++++++++++++ .../Bremsstrahlung/BremsstrahlungFunc.cpp | 66 +++ .../Bremsstrahlung/CMakeLists.txt | 8 + .../Bremsstrahlung/Make.package | 4 + .../Bremsstrahlung/PhotonCreationFunc.H | 212 ++++++++ .../Bremsstrahlung/PhotonCreationFunc.cpp | 35 ++ .../Collision/BinaryCollision/CMakeLists.txt | 1 + .../Coulomb/PairWiseCoulombCollisionFunc.H | 2 + .../Collision/BinaryCollision/DSMC/DSMCFunc.H | 3 + .../DSMC/SplitAndScatterFunc.H | 3 +- .../Collision/BinaryCollision/Make.package | 1 + .../NuclearFusion/NuclearFusionFunc.H | 3 + .../BinaryCollision/ParticleCreationFunc.H | 6 +- .../Particles/Collision/CollisionHandler.cpp | 8 + Source/Utils/ParticleUtils.H | 35 ++ 18 files changed, 880 insertions(+), 9 deletions(-) create mode 100644 Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H create mode 100644 Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp create mode 100644 Source/Particles/Collision/BinaryCollision/Bremsstrahlung/CMakeLists.txt create mode 100644 Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package create mode 100644 Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H create mode 100644 Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index b64c6d4b1fa..b4eaadeee76 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -11,6 +11,7 @@ #include "Particles/Collision/BinaryCollision/Coulomb/ComputeTemperature.H" #include "Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H" #include "Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H" +#include "Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H" #include "Particles/Collision/BinaryCollision/ParticleCreationFunc.H" #include "Particles/Collision/BinaryCollision/ShuffleFisherYates.H" #include "Particles/Collision/CollisionBase.H" @@ -363,6 +364,10 @@ public: amrex::Gpu::DeviceVector pair_reaction_weight(n_total_pairs); amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight = pair_reaction_weight.dataPtr(); + // Extra data needed when products are created + int const n_product_data = (binary_collision_functor.m_need_product_data ? n_total_pairs : 0); + amrex::Gpu::DeviceVector product_data(n_product_data); + amrex::ParticleReal* AMREX_RESTRICT p_product_data = product_data.dataPtr(); /* End of calculations only required when creating product particles */ @@ -458,7 +463,7 @@ public: // Call the function in order to perform collisions // If there are product species, mask, p_pair_indices_1/2, and - // p_pair_reaction_weight are filled here + // p_pair_reaction_weight and p_product_data are filled here binary_collision_functor( cell_start_1, cell_half_1, cell_half_1, cell_stop_1, @@ -467,7 +472,7 @@ public: n1, n1, T1, T1, q1, q1, m1, m1, dt, dV*volume_factor(i_cell), coll_idx, cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2, - p_pair_reaction_weight, engine); + p_pair_reaction_weight, p_product_data, engine); } ); @@ -481,7 +486,7 @@ public: products_mass, p_mask, products_np, copy_species1, copy_species2, p_pair_indices_1, p_pair_indices_2, - p_pair_reaction_weight); + p_pair_reaction_weight, p_product_data); for (int i = 0; i < n_product_species; i++) { @@ -589,6 +594,10 @@ public: amrex::Gpu::DeviceVector pair_reaction_weight(n_total_pairs); amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight = pair_reaction_weight.dataPtr(); + // Extra data needed when products are created + int const n_product_data = (binary_collision_functor.m_need_product_data ? n_total_pairs : 0); + amrex::Gpu::DeviceVector product_data(n_product_data); + amrex::ParticleReal* AMREX_RESTRICT p_product_data = product_data.dataPtr(); /* End of calculations only required when creating product particles */ @@ -717,7 +726,7 @@ public: // Call the function in order to perform collisions // If there are product species, p_mask, p_pair_indices_1/2, and - // p_pair_reaction_weight are filled here + // p_pair_reaction_weight and p_product_data are filled here binary_collision_functor( cell_start_1, cell_stop_1, cell_start_2, cell_stop_2, indices_1, indices_2, @@ -725,7 +734,7 @@ public: n1, n2, T1, T2, q1, q2, m1, m2, dt, dV*volume_factor(i_cell), coll_idx, cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2, - p_pair_reaction_weight, engine); + p_pair_reaction_weight, p_product_data, engine); } ); @@ -739,7 +748,7 @@ public: products_mass, p_mask, products_np, copy_species1, copy_species2, p_pair_indices_1, p_pair_indices_2, - p_pair_reaction_weight); + p_pair_reaction_weight, p_product_data); for (int i = 0; i < n_product_species; i++) { diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H index 2e24a93a35e..4341a6a496f 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.H @@ -19,6 +19,7 @@ enum struct CollisionType { DeuteriumTritiumToNeutronHeliumFusion, DeuteriumDeuteriumToNeutronHeliumFusion, DeuteriumHeliumToProtonHeliumFusion, ProtonBoronToAlphasFusion, + Bremsstrahlung, DSMC, PairwiseCoulomb, Undefined }; diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp index 05756583554..2740873ecba 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp @@ -31,6 +31,9 @@ namespace BinaryCollisionUtils{ const NuclearFusionType fusion_type = get_nuclear_fusion_type(collision_name, mypc); return nuclear_fusion_type_to_collision_type(fusion_type); } + else if (type == "bremsstrahlung") { + return CollisionType::Bremsstrahlung; + } else if (type == "dsmc") { return CollisionType::DSMC; } diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H new file mode 100644 index 00000000000..a21ef31e72d --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -0,0 +1,477 @@ +/* Copyright 2024 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_BREMSSTRAHLUNG_FUNC_H_ +#define WARPX_BREMSSTRAHLUNG_FUNC_H_ + +#include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H" +#include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/MultiParticleContainer.H" +#include "Particles/WarpXParticleContainer.H" +#include "Utils/Parser/ParserUtils.H" +#include "Utils/ParticleUtils.H" +#include "Utils/TextMsg.H" +#include "WarpX.H" + +#include +#include +#include +#include +#include +#include + +/** + * \brief This functor does binary Bremsstrahlung reactions in a single cell. + * Particles of the two interacting species are paired with each other and for each pair we compute + * the energy of the resulting photon if Bremsstrahlung happens. + * We fill p_mask with true if it happens so that the new photon corresponding to a given pair can be + * effectively created in the particle creation functor. + */ +class BremsstrahlungFunc +{ + // Define shortcuts for frequently-used type names + using ParticleType = WarpXParticleContainer::ParticleType; + using ParticleTileType = WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; + using index_type = ParticleBins::index_type; + using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + +public: + /* + * \brief Default constructor of the BremsstrahlungFunc class. + */ + BremsstrahlungFunc () = default; + + /** + * \brief Constructor of the BremsstrahlungFunc class + * + * @param[in] collision_name the name of the collision + * @param[in] mypc pointer to the MultiParticleContainer + * @param[in] isSameSpecies whether the two colliding species are the same (should always be false) + */ + BremsstrahlungFunc (std::string const& collision_name, MultiParticleContainer const * const mypc, + bool /*isSameSpecies*/); + + struct Executor { + /** + * \brief Executor of the BremsstrahlungFunc class. Produces Bremsstrahlung photons at the cell level. + * Note that this function does not yet create the resulting photons, but + * instead sets p_mask and calculates the weight and energy + * + * @param[in] I1s,I2s is the start index for I1,I2 (inclusive). + * @param[in] I1e,I2e is the stop index for I1,I2 (exclusive). + * @param[in] I1,I2 index arrays. They determine all elements that will be used. + * @param[in] soa_1,soa_2 contain the struct of array data of the two species + * @param[in] m1,m2 are masses. + * @param[in] dt is the time step length between two collision calls. + * @param[in] dV is the volume of the corresponding cell. + * @param[in] coll_idx is the collision index offset. + * @param[in] cell_start_pair is the start index of the pairs in that cell. + * @param[out] p_mask is a mask that will be set to true for all pairs. + * @param[out] p_pair_indices_1,p_pair_indices_2 arrays that store the indices of the + * particles of a given pair. Only p_pair_indices_1 will be used when creating photons. + * @param[out] p_pair_reaction_weight stores the weight of the resulting photons. + * @param[out] p_product_data stores the energy of the resulting photons. + * @param[in] engine the random engine. + */ + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void operator() ( + index_type const I1s, index_type const I1e, + index_type const I2s, index_type const I2e, + index_type const* AMREX_RESTRICT I1, + index_type const* AMREX_RESTRICT I2, + const SoaData_type& soa_1, const SoaData_type& soa_2, + GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, + amrex::ParticleReal const n1, amrex::ParticleReal const n2, + amrex::ParticleReal const /*T1*/, amrex::ParticleReal const /*T2*/, + amrex::ParticleReal const /*q1*/, amrex::ParticleReal const /*q2*/, + amrex::ParticleReal const m1, amrex::ParticleReal const /*m2*/, + amrex::Real const dt, amrex::Real const /*dV*/, index_type coll_idx, + index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask, + index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2, + amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + amrex::ParticleReal* AMREX_RESTRICT p_product_data, + amrex::RandomEngine const& engine) const + { + amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; + amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; + amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; + amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz]; + + amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; + amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux]; + amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy]; + amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; + + // Number of macroparticles of each species + const index_type NI1 = I1e - I1s; + const index_type NI2 = I2e - I2s; + + // Only loop over the number of particles in the first species (the electrons) + const index_type max_N = NI1; + const index_type min_N = NI2; + + index_type pair_index = cell_start_pair + coll_idx; + +#if (defined WARPX_DIM_RZ) + amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta]; + amrex::ParticleReal * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta]; +#endif + + index_type i1 = I1s + coll_idx; + index_type i2 = I2s + coll_idx; + // we will start from collision number = coll_idx and then add + // stride (smaller set size) until we do all collisions (larger set size) + // If there are more atoms than electrons, this loop executes once + // If there are more electrons than atoms, this loops over the electrons matched to the atom + for (index_type k = coll_idx; k < max_N; k += min_N) + { + +#if (defined WARPX_DIM_RZ) + /* In RZ geometry, macroparticles can collide with other macroparticles + * in the same *cylindrical* cell. For this reason, collisions between macroparticles + * are actually not local in space. In this case, the underlying assumption is that + * particles within the same cylindrical cell represent a cylindrically-symmetry + * momentum distribution function. Therefore, here, we temporarily rotate the + * momentum of one of the macroparticles in agreement with this cylindrical symmetry. + * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation; + * there is a corresponding assert statement at initialization.) */ + amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]]; + amrex::ParticleReal const u1xbuf = u1x[I1[i1]]; + u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta); + u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta); +#endif + + BremsstrahlungEvent( + u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], + u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], + m1, w1[ I1[i1] ], w2[ I2[i2] ], + n1, n2, dt, static_cast(pair_index), p_mask, + p_pair_reaction_weight, p_product_data, + engine); + +#if (defined WARPX_DIM_RZ) + amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]]; + u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta); + u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta); +#endif + + } + + p_pair_indices_1[pair_index] = I1[i1]; + p_pair_indices_2[pair_index] = I2[i2]; + i1 += min_N; + pair_index += min_N; + } + + static int constexpr nkoT1 = 17; + static int constexpr nKE = 11; + + /* \brief Calculate the cross section given the electron energy. + * The differential cross section is cut off below the plasma frequency since the + * low energy photons would all be immediately absorbed by the surrounding plasma. + * The normalized sigmaC is returned to allow calculating the photon energy. + * + * \param[in] E electron energy in Joules + * \param[in] n1 the electron number density + * \param[out] sigmaC the normalized differential cross section + * \result the cross section + */ + AMREX_GPU_HOST_DEVICE AMREX_INLINE + amrex::ParticleReal + CalculateCrossSection (amrex::ParticleReal KErel, + amrex::ParticleReal n1, + amrex::ParticleReal m1, + amrex::GpuArray & sigmaC, + int & i0_cut, + amrex::ParticleReal & koT1_cut) const + { + using namespace amrex::literals; + + amrex::ParticleReal const KErel_eV = KErel/PhysConst::q_e; + amrex::ParticleReal const wpe = PhysConst::q_e*std::sqrt(n1/(m1*PhysConst::ep0)); + amrex::ParticleReal const KEcut_eV = PhysConst::hbar*wpe/PhysConst::q_e; + + if (KErel_eV < KEcut_eV) { + // Ignore low energy electrons which would only produce low energy photons + return 0._prt; + } + + // find lo-index for koT1 cutoff (will typically be 1) + koT1_cut = std::max(KEcut_eV/KErel_eV, 1.0e-4); + i0_cut = 0; + for (int i=1; i < nkoT1; i++) { + if (m_koT1_grid[i] > koT1_cut) { + break; + } else { + i0_cut = i; + } + } + if (i0_cut == nkoT1 - 1) { + return 0.0_prt; + } + + // kdsigdk is linearly weighted in electron energy + amrex::ParticleReal w0, w1; + int index; + if (KErel_eV >= m_KEgrid_eV[nkoT1-1]) { + index = nkoT1 - 2; + w0 = 0._prt; + w1 = 1._prt; + } else { + // find index for electron energy interpolation + index = amrex::bisect(m_KEgrid_eV.data(), 0, nKE, KErel_eV); + + // compute interpolation weights for k*dsigma/dk + w0 = (m_KEgrid_eV[index+1] - KErel_eV)/(m_KEgrid_eV[index+1] - m_KEgrid_eV[index]); + w1 = 1.0 - w0; + } + + // kdsigdk at the cutoff + amrex::ParticleReal const w00 = (m_koT1_grid[i0_cut+1] - koT1_cut)/(m_koT1_grid[i0_cut+1] - m_koT1_grid[i0_cut]); + amrex::ParticleReal const w01 = 1.0 - w00; + amrex::ParticleReal const kdsigdk_cut = ( w00*(w0*m_kdsigdk[index][i0_cut] + w1*m_kdsigdk[index+1][i0_cut]) + + w01*(w0*m_kdsigdk[index][i0_cut+1] + w1*m_kdsigdk[index+1][i0_cut+1])); + + // create cumulative distribution using trapezoidal rule + amrex::ParticleReal m_koT1_grid_im1 = koT1_cut; + amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; + for (int i=i0_cut+1; i < nkoT1; i++) { + amrex::ParticleReal const this_dk = (m_koT1_grid[i] - m_koT1_grid_im1); + amrex::ParticleReal const this_k = (m_koT1_grid[i] + m_koT1_grid_im1)*0.5_prt; + amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + w1*m_kdsigdk[index+1][i]; + // using centered k here mitigates divergece effect for k->0 + amrex::ParticleReal const this_dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5/this_k; + sigmaC[i] = sigmaC[i-1] + this_dsigdk*this_dk; + m_koT1_grid_im1 = m_koT1_grid[i]; + kdsigdk_im1 = kdsigdk_i; + } + amrex::ParticleReal const sigma_total = sigmaC[nkoT1-1]; // total cross section [m^2] + + for (int i=i0_cut+1; i < nkoT1; i++) { + sigmaC[i] /= sigma_total; + } + + return sigma_total; + } + + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void BremsstrahlungEvent (amrex::ParticleReal & u1x, amrex::ParticleReal & u1y, + amrex::ParticleReal & u1z, amrex::ParticleReal const& u2x, + amrex::ParticleReal const& u2y, amrex::ParticleReal const& u2z, + amrex::ParticleReal m1, + amrex::ParticleReal w1, amrex::ParticleReal /*w2*/, + amrex::ParticleReal n1, amrex::ParticleReal n2, + amrex::Real dt, int pair_index, + index_type* AMREX_RESTRICT p_mask, + amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + amrex::ParticleReal* AMREX_RESTRICT p_product_data, + amrex::RandomEngine const& engine) const + { + using namespace amrex::literals; + + constexpr auto c2 = PhysConst::c * PhysConst::c; + amrex::ParticleReal const m_e_eV = m1*PhysConst::c*PhysConst::c/PhysConst::q_e; // 0.511e6 + + // compute gamma for particles 1 (electron) and 2 (ion) + amrex::ParticleReal gbp1sq = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; + amrex::ParticleReal gbp2sq = (u2x*u2x + u2y*u2y + u2z*u2z)/c2; + amrex::ParticleReal gammap1 = std::sqrt(1.0_prt + gbp1sq); + amrex::ParticleReal gammap2 = std::sqrt(1.0_prt + gbp2sq); + + // transform electron to frame where ion is at rest + amrex::ParticleReal u1x_rel = u1x; + amrex::ParticleReal u1y_rel = u1y; + amrex::ParticleReal u1z_rel = u1z; + ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, u2x, u2y, u2z); + + // compute electron KE in frame where ion is at rest + amrex::ParticleReal gbp1sq_rel = (u1x_rel*u1x_rel + u1y_rel*u1y_rel + u1z_rel*u1z_rel)/c2; + amrex::ParticleReal gammap1_rel = std::sqrt(1.0_prt + gbp1sq_rel); + amrex::ParticleReal KE_eV = m_e_eV*gbp1sq_rel/(1.0_prt + gammap1_rel); + + amrex::GpuArray sigmaC; + sigmaC.fill(0._prt); + int i0_cut; + amrex::ParticleReal koT1_cut; + amrex::ParticleReal const sigma_total = CalculateCrossSection(KE_eV, n1, m1, sigmaC, i0_cut, koT1_cut); + + if (sigma_total == 0._prt) { return; } + + // determine if the pair collide and then do the collision + amrex::ParticleReal fmulti = m_multiplier; // >= 1.0 + amrex::ParticleReal const gamma_factor = gammap1_rel/(gammap1*gammap2); + amrex::ParticleReal const beta1 = std::sqrt(gbp1sq_rel)/gammap1_rel; // magnitude electron beta + amrex::ParticleReal arg = fmulti*beta1*PhysConst::c*sigma_total*n2*dt*gamma_factor; + if (arg > 1.0_prt) { + amrex::Print() << "Notice: arg = " << arg << " in interSpeciesBremsstrahlung" << "\n"; + amrex::Print() << " beta1 = " << beta1 << "\n"; + amrex::Print() << " n2 = " << n2 << "\n"; + amrex::Print() << " sigma_total = " << sigma_total << "\n"; + arg = 1.0_prt; + fmulti = fmulti/arg; + if (fmulti < 1.0_prt) { fmulti = 1.0_prt; } + } + //q12 = 1.0 - exp(-arg); + amrex::ParticleReal q12 = arg; + + // Get a random number + amrex::ParticleReal const random_number = amrex::Random(engine); + p_mask[pair_index] = (random_number < q12); + if (random_number < q12) { + + // compute weight for photon + amrex::ParticleReal const w_photon = w1/fmulti; + + // compute lab-frame electron kinetic energy before emission + amrex::ParticleReal KE_before = m_e_eV*gbp1sq/(1.0_prt + gammap1); + + // get energy of photon (hbar*omega) + amrex::ParticleReal const random_number2 = amrex::Random(engine); + int const index = amrex::bisect(sigmaC.data(), i0_cut, nkoT1, random_number2); + amrex::ParticleReal const w0 = (random_number2 - sigmaC[index])/(sigmaC[index+1] - sigmaC[index]); + amrex::ParticleReal const Epi = (index == i0_cut ? koT1_cut : m_koT1_grid[index]); + amrex::ParticleReal const EphotonoT1 = (1._prt - w0)*Epi + w0*m_koT1_grid[index+1]; + amrex::ParticleReal const Ephoton_eV = EphotonoT1*KE_eV; + + // adjust electron gamma and beta + amrex::ParticleReal const KE_prime_eV = KE_eV - Ephoton_eV*w_photon/w1; + amrex::ParticleReal const gamma_prime = 1.0_prt + KE_prime_eV/m_e_eV; + + // rescale electron beta in frame where ion is at rest + amrex::ParticleReal const u_rel_ratio = std::sqrt((gamma_prime + 1.0_prt)*KE_prime_eV/m_e_eV/gbp1sq_rel); + u1x_rel *= u_rel_ratio; + u1y_rel *= u_rel_ratio; + u1z_rel *= u_rel_ratio; + + // Lorentz transform electron back to lab frame + ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, -u2x, -u2y, -u2z); + u1x = u1x_rel; + u1y = u1y_rel; + u1z = u1z_rel; + + // compute electron kinetic energy after emission + gbp1sq = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; + gammap1 = std::sqrt(1.0_prt + gbp1sq); + amrex::ParticleReal const KE_after = m_e_eV*gbp1sq/(1.0_prt + gammap1); + + // update energy lost to bremsstrahlung + amrex::ParticleReal const deltaE_eV = KE_after - KE_before; + /* m_deltaE_bremsstrahlung -= deltaE_eV*w1*PhysConst::q_e; // Joules */ + + // compute photon energy in lab frame + amrex::ParticleReal const Ephoton_lab = deltaE_eV*w1/w_photon*PhysConst::q_e; // Joules + + p_pair_reaction_weight[pair_index] = w_photon; + p_product_data[pair_index] = Ephoton_lab; + + } + + } + + bool m_computeSpeciesDensities = true; + bool m_computeSpeciesTemperatures = false; + bool m_need_product_data = true; + amrex::ParticleReal m_multiplier; + + // relative photon energy grid ( koT1 = Ephoton_eV/KErel_eV ) + amrex::GpuArray m_koT1_grid = {0., 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.99, 1.00}; + + // energy grid for electrons in eV + amrex::GpuArray m_KEgrid_eV = {1.0e3, 2.0e3, 5.0e3, 1.0e4, 2.0e4, 5.0e4, 1.0e5, 2.0e5, 5.0e5, 1.0e6, 2.0e6}; + + // differential total cross section + amrex::GpuArray, nKE> m_kdsigdk; + + }; + + [[nodiscard]] Executor const& executor () const { return m_exe; } + +private: + + Executor m_exe; + + // Cross section data + + void UploadCrossSection (int Z); + + // scaled energy-weighted differential total cross section for Bremsstrahlung (kdsigdk_n + Z*kdsigdk_e) + // for e + H from Table I of Seltzer and Berger, ATOMIC DATA AND NUCLEAR DATA TABLES 35,345-418 (1985). + // The values here are actualy beta1^2/Z^2*k*dsig/dk in units of mBarns = 1e-31 m^2 where Z is the + // atomic number and beta1 = sqrt(1.0 - 1/gamma1^2) with gamma1 = 1 + T1/(me*c^2). These vectors are + // converted to k*dsig/dk in units of m^2 on initialization. + + std::map>> m_kdsigdk_map = { + + // atomic number Z = 1 + {1, { + {7.853, 7.849, 7.833, 7.800, 7.746, 7.446, 7.040, 6.586, 6.124, 5.664, 5.230, 4.841, 4.521, 4.400, 4.372, 4.362, 4.360}, // 1.0 keV + {8.805, 8.817, 8.801, 8.738, 8.638, 8.059, 7.377, 6.699, 6.052, 5.431, 4.839, 4.263, 3.682, 3.437, 3.374, 3.326, 3.320}, // 2.0 keV + {10.32, 10.25, 10.15, 9.976, 9.753, 8.703, 7.661, 6.728, 5.899, 5.148, 4.424, 3.697, 2.897, 2.446, 2.286, 2.161, 2.104}, // 5.0 keV + {11.55, 11.40, 11.19, 10.89, 10.55, 9.087, 7.795, 6.711, 5.776, 4.952, 4.172, 3.387, 2.509, 1.982, 1.762, 1.548, 1.439}, // 10 keV + {13.07, 12.61, 12.13, 11.62, 11.11, 9.370, 7.915, 6.674, 5.617, 4.734, 3.933, 3.141, 2.239, 1.684, 1.424, 1.129, 0.967}, // 20 keV + {14.90, 14.17, 13.42, 12.65, 11.92, 9.620, 7.858, 6.426, 5.295, 4.370, 3.575, 2.777, 1.940, 1.401, 1.112, 0.754, 0.551}, // 50 keV + {16.94, 15.72, 14.55, 13.47, 12.51, 9.686, 7.689, 6.107, 4.937, 3.951, 3.179, 2.438, 1.670, 1.180, 0.903, 0.548, 0.343}, // 100 keV + {20.31, 18.17, 16.25, 14.71, 13.46, 9.837, 7.509, 5.808, 4.595, 3.549, 2.694, 2.022, 1.362, 0.938, 0.697, 0.380, 0.197}, // 200 keV + {27.90, 23.39, 19.68, 17.28, 15.68, 10.84, 7.831, 5.945, 4.555, 3.400, 2.386, 1.634, 1.045, 0.706, 0.505, 0.244, 0.094}, // 500 keV + {35.57, 28.83, 23.45, 20.32, 18.41, 12.21, 8.960, 6.849, 5.241, 3.927, 2.763, 1.732, 1.038, 0.688, 0.483, 0.218, 0.065}, // 1.0 MeV + {43.42, 34.55, 27.56, 23.66, 21.45, 14.61, 11.03, 8.605, 6.744, 5.195, 3.809, 2.472, 1.306, 0.846, 0.587, 0.250, 0.057} // 2.0 MeV + }}, + + // atomic number Z = 2 + {2, { + {7.167, 7.192, 7.206, 7.201, 7.181, 7.001, 6.726, 6.409, 6.077, 5.740, 5.422, 5.150, 4.939, 4.858, 4.837, 4.825, 4.821}, // 1.0 keV + {8.232, 8.239, 8.229, 8.187, 8.120, 7.713, 7.200, 6.666, 6.147, 5.646, 5.173, 4.729, 4.313, 4.145, 4.099, 4.064, 4.053}, // 2.0 keV + {9.678, 9.640, 9.570, 9.444, 9.276, 8.439, 7.568, 6.770, 6.055, 5.397, 4.770, 4.135, 3.510, 3.180, 3.070, 2.986, 2.951}, // 5.0 keV + {10.81, 10.71, 10.56, 10.33, 10.06, 8.834, 7.706, 6.747, 5.914, 5.166, 4.461, 3.762, 3.009, 2.590, 2.430, 2.283, 2.213}, // 10 keV + {12.18, 11.80, 11.40, 10.98, 10.55, 9.062, 7.784, 6.678, 5.721, 4.897, 4.148, 3.420, 2.605, 2.131, 1.927, 1.711, 1.600}, // 20 keV + {13.49, 12.92, 12.33, 11.71, 11.11, 9.139, 7.592, 6.336, 5.325, 4.473, 3.707, 2.950, 2.148, 1.657, 1.413, 1.130, 0.974}, // 50 keV + {14.54, 13.71, 12.88, 12.05, 11.28, 8.928, 7.241, 5.917, 4.892, 4.017, 3.265, 2.544, 1.797, 1.334, 1.093, 0.792, 0.624}, // 100 keV + {16.12, 14.79, 13.54, 12.44, 11.49, 8.747, 6.868, 5.470, 4.418, 3.513, 2.751, 2.092, 1.433, 1.023, 0.802, 0.525, 0.365}, // 200 keV + {19.94, 17.44, 16.27, 13.69, 12.50, 9.013, 6.783, 5.288, 4.140, 3.184, 2.345, 1.667, 1.080, 0.745, 0.556, 0.315, 0.178}, // 500 keV + {24.04, 20.39, 17.38, 15.42, 14.06, 9.865, 7.470, 5.820, 4.535, 3.479, 2.549, 1.714, 1.060, 0.713, 0.517, 0.266, 0.123}, // 1.0 MeV + {28.11, 23.46, 19.71, 17.43, 15.98, 11.46, 8.864, 7.025, 5.582, 4.375, 3.303, 2.268, 1.320, 0.866, 0.613, 0.292, 0.108} // 2.0 MeV + }}, + + // atomic number Z = 5 + {5, { + {5.670, 5.717, 5.760, 5.793, 5.820, 5.871, 5.856, 5.797, 5.713, 5.610, 5.502, 5.411, 5.350, 5.323, 5.311, 5.299, 5.293}, // 1.0 keV + {6.814, 6.863, 6.903, 6.925, 6.932, 6.863, 6.698, 6.484, 6.250, 6.010, 5.789, 5.596, 5.438, 5.366, 5.340, 5.318, 5.306}, // 2.0 keV + {8.325, 8.350, 8.360, 8.324, 8.265, 7.870, 7.382, 6.895, 6.441, 6.020, 5.635, 5.285, 4.987, 4.880, 4.847, 4.814, 4.799}, // 5.0 keV + {9.457, 9.432, 9.377, 9.272, 9.128, 8.384, 7.614, 6.920, 6.310, 5.759, 5.259, 4.795, 4.372, 4.201, 4.150, 4.111, 4.097}, // 10 keV + {10.56, 10.44, 10.29, 10.07, 9.819, 8.684, 7.651, 6.784, 6.042, 5.381, 4.785, 4.222, 3.676, 3.414, 3.332, 3.275, 3.252}, // 20 keV + {12.00, 11.56, 11.11, 10.67, 10.23, 8.718, 7.455, 6.396, 5.515, 4.765, 4.097, 3.463, 2.785, 2.431, 2.302, 2.185, 2.136}, // 50 keV + {12.69, 12.07, 11.45, 10.86, 10.30, 8.452, 7.031, 5.897, 4.996, 4.223, 3.526, 2.867, 2.182, 1.807, 1.650, 1.490, 1.410}, // 100 keV + {13.49, 12.57, 11.70, 10.91, 10.22, 8.039, 6.503, 5.341, 4.421, 3.625, 2.914, 2.286, 1.651, 1.287, 1.121, 0.940, 0.844}, // 200 keV + {15.27, 13.74, 12.38, 11.33, 10.49, 7.918, 6.188, 4.945, 3.953, 3.121, 2.392, 1.765, 1.185, 0.865, 0.706, 0.521, 0.419}, // 500 keV + {17.23, 15.27, 13.58, 12.36, 11.41, 8.400, 6.550, 5.235, 4.166, 3.260, 2.464, 1.749, 1.131, 0.788, 0.617, 0.409, 0.293}, // 1.0 MeV + {19.34, 16.89, 14.86, 13.50, 12.53, 9.504, 7.545, 6.077, 4.899, 3.908, 3.027, 2.183, 1.368, 0.924, 0.693, 0.414, 0.256} // 2.0 MeV + }}, + + // atomic number Z = 6 + {6, { + {5.336, 5.384, 5.427, 5.464, 5.494, 5.570, 5.585, 5.557, 5.501, 5.425, 5.337, 5.259, 5.206, 5.185, 5.175, 5.164, 5.158}, // 1.0 keV + {6.498, 6.553, 6.600, 6.630, 6.648, 6.628, 6.518, 6.359, 6.174, 5.976, 5.790, 5.619, 5.470, 5.400, 5.375, 5.355, 5.346}, // 2.0 keV + {8.028, 8.065, 8.084, 8.070, 8.031, 7.721, 7.315, 6.897, 6.502, 6.135, 5.801, 5.503, 5.251, 5.160, 5.129, 5.096, 5.080}, // 5.0 keV + {9.168, 9.159, 9.123, 9.041, 8.922, 8.281, 7.593, 6.964, 6.411, 5.913, 5.467, 5.064, 4.712, 4.577, 4.536, 4.500, 4.485}, // 10 keV + {10.27, 10.18, 10.05, 9.865, 9.638, 8.606, 7.644, 6.832, 6.140, 5.522, 4.973, 4.465, 3.994, 3.779, 3.715, 3.671, 3.654}, // 20 keV + {11.72, 11.31, 10.90, 10.49, 10.08, 8.651, 7.450, 6.436, 5.587, 4.862, 4.221, 3.622, 2.997, 2.684, 2.578, 2.489, 2.455}, // 50 keV + {12.38, 11.81, 11.24, 10.68, 10.16, 8.400, 7.026, 5.924, 5.045, 4.293, 3.611, 2.971, 2.314, 1.967, 1.830, 1.701, 1.638}, // 100 keV + {13.10, 12.26, 11.45, 10.72, 10.06, 7.971, 6.481, 5.348, 4.447, 3.668, 2.970, 2.351, 1.726, 1.378, 1.225, 1.068, 0.987}, // 200 keV + {14.65, 13.27, 12.04, 11.06, 10.26, 7.806, 6.132, 4.921, 3.953, 3.133, 2.417, 1.797, 1.222, 0.904, 0.756, 0.586, 0.494}, // 500 keV + {16.39, 14.64, 13.12, 11.99, 11.09, 8.243, 6.460, 5.179, 4.135, 3.248, 2.467, 1.769, 1.153, 0.814, 0.650, 0.453, 0.345}, // 1.0 MeV + {18.14, 16.07, 14.31, 13.09, 12.18, 9.254, 7.367, 5.974, 4.837, 3.869, 3.000, 2.184, 1.385, 0.942, 0.720, 0.453, 0.301} // 2.0 MeV + }} + + }; + +}; + +#endif // WARPX_BREMSSTRAHLUNG_FUNC_H_ + + diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp new file mode 100644 index 00000000000..5945127d013 --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -0,0 +1,66 @@ +/* Copyright 2024 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "BremsstrahlungFunc.H" + +#include +#include + +using namespace amrex::literals; + +BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, MultiParticleContainer const * const mypc, + bool /*isSameSpecies*/) +{ + const amrex::ParmParse pp_collision_name(collision_name); + + // Read in the number of electrons on the target + int Z; + pp_collision_name.get("Z", Z); + + std::string product_species_name; + pp_collision_name.get("product_species", product_species_name); + auto& product_species = mypc->GetParticleContainerFromName(product_species_name); + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(product_species.AmIA(), + "The product species must be photons"); + + amrex::ParticleReal multiplier = 1._prt; + pp_collision_name.get("multiplier", multiplier); + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(multiplier >= 1., + "The multiplier must be greater than or equal to one"); + m_exe.m_multiplier = multiplier; + + // Fill in the m_kdsigdk array + UploadCrossSection(Z); + +} + +void +BremsstrahlungFunc::UploadCrossSection (int Z) +{ + + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(m_kdsigdk_map.count(Z) == 1, + "Bremsstrahlung cross section not available for Z = " + std::to_string(Z) + "!"); + + constexpr auto m_e_eV = PhysConst::m_e*PhysConst::c*PhysConst::c/PhysConst::q_e; // 0.511e6 + + std::vector> & kdsigdk = m_kdsigdk_map.at(2); + + // Convert Seltzer and Berger energy-weighted differential cross section to units of [m^2] + for (int iee=0; iee < m_exe.nKE; iee++) { + amrex::ParticleReal const E = m_exe.m_KEgrid_eV[iee]/m_e_eV; + amrex::ParticleReal const gamma = 1.0_rt + E; + /* betaSq = 1.0 - 1.0/gamma/gamma */ + amrex::ParticleReal const betaSq = (E*E + 2._rt*E)/gamma/gamma; + amrex::ParticleReal const scale_factor = 1.0e-31_rt*Z*Z/betaSq; + for (int iep=0; iep < m_exe.nkoT1; iep++) { + m_exe.m_kdsigdk[iee][iep] = kdsigdk[iee][iep]*scale_factor; + } + } + +} + diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/CMakeLists.txt b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/CMakeLists.txt new file mode 100644 index 00000000000..8e0aa37b5a4 --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/CMakeLists.txt @@ -0,0 +1,8 @@ +foreach(D IN LISTS WarpX_DIMS) + warpx_set_suffix_dims(SD ${D}) + target_sources(lib_${SD} + PRIVATE + BremsstrahlungFunc.cpp + PhotonCreationFunc.cpp + ) +endforeach() diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package new file mode 100644 index 00000000000..7bc4c15591c --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package @@ -0,0 +1,4 @@ +CEXE_sources += BremsstrahlungFunc.cpp +CEXE_sources += PhotonCreationFunc.cpp + +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision/Bremmstrahlung diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H new file mode 100644 index 00000000000..ac9c213b104 --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H @@ -0,0 +1,212 @@ +/* Copyright 2024 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef WARPX_PHOTON_CREATION_FUNC_H_ +#define WARPX_PHOTON_CREATION_FUNC_H_ + +#include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H" + +#include "Particles/ParticleCreation/SmartCopy.H" +#include "Particles/MultiParticleContainer.H" +#include "Particles/WarpXParticleContainer.H" + +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * \brief This functor creates photons produced from a binary collision and sets their initial + * properties (position, momentum, weight). + */ +class PhotonCreationFunc +{ + // Define shortcuts for frequently-used type names + using ParticleType = typename WarpXParticleContainer::ParticleType; + using ParticleTileType = typename WarpXParticleContainer::ParticleTileType; + using ParticleTileDataType = typename ParticleTileType::ParticleTileDataType; + using ParticleBins = amrex::DenseBins; + using index_type = typename ParticleBins::index_type; + using SoaData_type = typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + +public: + /** + * \brief Default constructor of the PhotonCreationFunc class. + */ + PhotonCreationFunc () = default; + + /** + * \brief Constructor of the PhotonCreationFunc class + * + * @param[in] collision_name the name of the collision + * @param[in] mypc pointer to the MultiParticleContainer + */ + PhotonCreationFunc (const std::string& collision_name, MultiParticleContainer const * mypc); + + /** + * \brief operator() of the PhotonCreationFunc class. It creates new photon particles from binary + * collisions. + * One product particle is created at the position of the first parent particle that collided. + * For example, with the Bremsstrahlung collision between an electron and an ion, + * the new photon is created at the position of the electron. + * This function also sets the initial weight of the new photon to be the same as the weight + * of the first parent. + * This function sets the initial momentum of the product particles, by calling a + * function specific to the considered binary collision. + * Finally, this function subtracts the energy of the new photon from the first parent, i.e. the electron. + * + * @param[in] n_total_pairs how many binary collisions have been performed in this tile + * @param[in, out] ptile1,ptile2 the particle tiles of the two colliding species + * @param[out] tile_products array containing tile data of the product particles. + * @param[in] m1 mass of the first colliding particle species + * @param[in] m2 mass of the second colliding particle species + * @param[in] products_mass array storing the mass of product particles + * @param[in] p_mask a mask that is 1 if binary collision has resulted in particle creation + * event, 0 otherwise. + * @param[in] products_np array storing the number of existing product particles in that tile + * @param[in] copy_species1 array of functors used to copy data from the first colliding + * particle species to the product particles and to default initialize product particle + * quantities. + * @param[in] copy_species2 array of functors used to copy data from the second colliding + * particle species to the product particles and to default initialize product particle + * quantities. + * @param[in] p_pair_indices_1 array that stores the indices of the particles of the first + * colliding species that were used in the binary collisions (i.e. particle with index + * p_pair_indices_1[i] took part in collision i) + * @param[in] p_pair_indices_2 array that stores the indices of the particles of the second + * colliding species that were used in the binary collisions (i.e. particle with index + * p_pair_indices_2[i] took part in collision i) + * @param[in] p_pair_reaction_weight array that stores the photon weight + * @param[in] p_product_data array that stores the photon energy + */ + AMREX_INLINE + amrex::Vector operator() ( + const index_type& n_total_pairs, + ParticleTileType& ptile1, ParticleTileType& /*ptile2*/, + const amrex::Vector& pc_products, + ParticleTileType** AMREX_RESTRICT tile_products, + const amrex::ParticleReal& /*m1*/, const amrex::ParticleReal& /*m2*/, + const amrex::Vector& /*products_mass*/, + const index_type* AMREX_RESTRICT p_mask, + const amrex::Vector& products_np, + const SmartCopy* AMREX_RESTRICT copy_species1, + const SmartCopy* AMREX_RESTRICT /*copy_species2*/, + const index_type* AMREX_RESTRICT p_pair_indices_1, + const index_type* AMREX_RESTRICT /*p_pair_indices_2*/, + const amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + const amrex::ParticleReal* AMREX_RESTRICT p_product_data + ) const + { + using namespace amrex::literals; + + if (n_total_pairs == 0) { return amrex::Vector(m_num_product_species, 0); } + + // Compute offset array and allocate memory for the produced species + amrex::Gpu::DeviceVector offsets(n_total_pairs); + const auto total = amrex::Scan::ExclusiveSum(n_total_pairs, p_mask, offsets.data()); + const index_type* AMREX_RESTRICT p_offsets = offsets.dataPtr(); + + amrex::Vector num_added_vec(m_num_product_species); + for (int i = 0; i < m_num_product_species; i++) + { + const index_type num_added = total; + num_added_vec[i] = static_cast(num_added); + tile_products[i]->resize(products_np[i] + num_added); + } + + const auto soa_1 = ptile1.getParticleTileData(); + + // Create necessary GPU vectors, that will be used in the kernel below + amrex::Vector soa_products; + for (int i = 0; i < m_num_product_species; i++) + { + soa_products.push_back(tile_products[i]->getParticleTileData()); + } +#ifdef AMREX_USE_GPU + amrex::Gpu::DeviceVector device_soa_products(m_num_product_species); + amrex::Gpu::DeviceVector device_products_np(m_num_product_species); + + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, soa_products.begin(), + soa_products.end(), + device_soa_products.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, products_np.begin(), + products_np.end(), + device_products_np.begin()); + + amrex::Gpu::streamSynchronize(); + SoaData_type* AMREX_RESTRICT soa_products_data = device_soa_products.data(); + const index_type* AMREX_RESTRICT products_np_data = device_products_np.data(); +#else + SoaData_type* AMREX_RESTRICT soa_products_data = soa_products.data(); + const index_type* AMREX_RESTRICT products_np_data = products_np.data(); +#endif + + amrex::ParallelForRNG(n_total_pairs, + [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept + { + if (p_mask[i]) + { + const auto product_index = products_np_data[0] + p_offsets[i]; + + // Create product particle at position of particle 1 + copy_species1[0](soa_products_data[0], soa_1, static_cast(p_pair_indices_1[i]), + static_cast(product_index), engine); + + // Set the weight of the new particle to p_pair_reaction_weight[i] + soa_products_data[0].m_rdata[PIdx::w][product_index] = p_pair_reaction_weight[i]; + + auto& ux1 = soa_products_data[0].m_rdata[PIdx::ux][product_index]; + auto& uy1 = soa_products_data[0].m_rdata[PIdx::uy][product_index]; + auto& uz1 = soa_products_data[0].m_rdata[PIdx::uz][product_index]; + + // Normalize out the electron velocity and multiply by the photon momentum, E/c + auto u1 = std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1); + ux1 *= (p_product_data[i]/PhysConst::c)/u1; + uy1 *= (p_product_data[i]/PhysConst::c)/u1; + uz1 *= (p_product_data[i]/PhysConst::c)/u1; + } + }); + + // Initialize the user runtime components + for (int i = 0; i < m_num_product_species; i++) + { + const auto start_index = int(products_np[i]); + const auto stop_index = int(products_np[i] + num_added_vec[i]); + ParticleCreation::DefaultInitializeRuntimeAttributes(*tile_products[i], + 0, 0, + pc_products[i]->getUserRealAttribs(), pc_products[i]->getUserIntAttribs(), + pc_products[i]->getParticleComps(), pc_products[i]->getParticleiComps(), + pc_products[i]->getUserRealAttribParser(), + pc_products[i]->getUserIntAttribParser(), +#ifdef WARPX_QED + false, // do not initialize QED quantities, since they were initialized + // when calling the SmartCopy functors + pc_products[i]->get_breit_wheeler_engine_ptr(), + pc_products[i]->get_quantum_sync_engine_ptr(), +#endif + pc_products[i]->getIonizationInitialLevel(), + start_index, stop_index); + } + + amrex::Gpu::synchronize(); + + return num_added_vec; + } + +private: + // How many different type of species the collision produces + int m_num_product_species; + + CollisionType m_collision_type; + +}; + +#endif // WARPX_PHOTON_CREATION_FUNC_H_ diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp new file mode 100644 index 00000000000..120657a8382 --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp @@ -0,0 +1,35 @@ +/* Copyright 2024 David Grote + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#include "PhotonCreationFunc.H" + +#include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H" +#include "Particles/MultiParticleContainer.H" +#include "Utils/TextMsg.H" + +#include +#include +#include + +#include + +PhotonCreationFunc::PhotonCreationFunc (const std::string& collision_name, + MultiParticleContainer const * const mypc): + m_collision_type{BinaryCollisionUtils::get_collision_type(collision_name, mypc)} +{ + + if (m_collision_type == CollisionType::Bremsstrahlung) + { + + // Only photons created and only one photon per collision + m_num_product_species = 1; + + } else { + WARPX_ABORT_WITH_MESSAGE("Unknown collision type in PhotonCreationFunc"); + } + +} diff --git a/Source/Particles/Collision/BinaryCollision/CMakeLists.txt b/Source/Particles/Collision/BinaryCollision/CMakeLists.txt index 60933c83a21..6e5892cf829 100644 --- a/Source/Particles/Collision/BinaryCollision/CMakeLists.txt +++ b/Source/Particles/Collision/BinaryCollision/CMakeLists.txt @@ -8,3 +8,4 @@ foreach(D IN LISTS WarpX_DIMS) endforeach() add_subdirectory(DSMC) +add_subdirectory(Bremsstrahlung) diff --git a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H index 3322c19ed2d..42806a95f83 100644 --- a/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H @@ -98,6 +98,7 @@ public: index_type const /*cell_start_pair*/, index_type* /*p_mask*/, index_type* /*p_pair_indices_1*/, index_type* /*p_pair_indices_2*/, amrex::ParticleReal* /*p_pair_reaction_weight*/, + amrex::ParticleReal* /*p_product_data*/, amrex::RandomEngine const& engine) const { using namespace amrex::literals; @@ -112,6 +113,7 @@ public: amrex::ParticleReal m_CoulombLog; bool m_computeSpeciesDensities = true; bool m_computeSpeciesTemperatures = false; + bool m_need_product_data = false; bool m_isSameSpecies; }; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H index 5a3c925e9bd..40c5bae9506 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H @@ -86,6 +86,7 @@ public: * @param[out] p_pair_reaction_weight stores the weight of the product particles. It is only * needed here to store information that will be used later on when actually creating the * product particles. + * @param[out] p_product_data stores addition data for the product particles. It is not used here. * @param[in] engine the random engine. */ AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -104,6 +105,7 @@ public: index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask, index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2, amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + amrex::ParticleReal* AMREX_RESTRICT /*p_product_data*/, amrex::RandomEngine const& engine) const { amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; @@ -201,6 +203,7 @@ public: int m_process_count; bool m_computeSpeciesDensities = false; bool m_computeSpeciesTemperatures = false; + bool m_need_product_data = false; bool m_isSameSpecies = false; ScatteringProcess::Executor* m_scattering_processes_data; }; diff --git a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H index 473199a6b21..8aa528f6925 100644 --- a/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H +++ b/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H @@ -63,7 +63,8 @@ public: const SmartCopy* AMREX_RESTRICT copy_species2, const index_type* AMREX_RESTRICT p_pair_indices_1, const index_type* AMREX_RESTRICT p_pair_indices_2, - const amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight ) const + const amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + const amrex::ParticleReal* AMREX_RESTRICT /*p_product_data*/ ) const { using namespace amrex::literals; diff --git a/Source/Particles/Collision/BinaryCollision/Make.package b/Source/Particles/Collision/BinaryCollision/Make.package index 393e6270345..dde627b096f 100644 --- a/Source/Particles/Collision/BinaryCollision/Make.package +++ b/Source/Particles/Collision/BinaryCollision/Make.package @@ -2,5 +2,6 @@ CEXE_sources += BinaryCollisionUtils.cpp CEXE_sources += ParticleCreationFunc.cpp include $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision/DSMC/Make.package +include $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H index 14b1422db27..716b36f24b3 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -121,6 +121,7 @@ public: * @param[out] p_pair_reaction_weight stores the weight of the product particles. It is only * needed here to store information that will be used later on when actually creating the * product particles. + * @param[out] p_product_data stores addition data for the product particles. It is not used here. * @param[in] engine the random engine. */ AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -139,6 +140,7 @@ public: index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask, index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2, amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + amrex::ParticleReal* AMREX_RESTRICT /*p_product_data*/, amrex::RandomEngine const& engine) const { amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; @@ -239,6 +241,7 @@ public: NuclearFusionType m_fusion_type; bool m_computeSpeciesDensities = false; bool m_computeSpeciesTemperatures = false; + bool m_need_product_data = false; bool m_isSameSpecies; }; diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H index e4772aab7c9..c4df63139be 100644 --- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H @@ -104,7 +104,8 @@ public: const SmartCopy* AMREX_RESTRICT copy_species2, const index_type* AMREX_RESTRICT p_pair_indices_1, const index_type* AMREX_RESTRICT p_pair_indices_2, - const amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight + const amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + const amrex::ParticleReal* AMREX_RESTRICT /*p_product_data*/ ) const { using namespace amrex::literals; @@ -295,7 +296,8 @@ public: const index_type* /*p_mask*/, const amrex::Vector& /*products_np*/, const SmartCopy* /*copy_species1*/, const SmartCopy* /*copy_species2*/, const index_type* /*p_pair_indices_1*/, const index_type* /*p_pair_indices_2*/, - const amrex::ParticleReal* /*p_pair_reaction_weight*/ + const amrex::ParticleReal* /*p_pair_reaction_weight*/, + const amrex::ParticleReal* AMREX_RESTRICT /*p_product_data*/ ) const { return {}; diff --git a/Source/Particles/Collision/CollisionHandler.cpp b/Source/Particles/Collision/CollisionHandler.cpp index d723a1017d8..f7bee21e918 100644 --- a/Source/Particles/Collision/CollisionHandler.cpp +++ b/Source/Particles/Collision/CollisionHandler.cpp @@ -13,6 +13,8 @@ #include "Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H" #include "Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H" #include "Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H" +#include "Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H" +#include "Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H" #include "Particles/Collision/BinaryCollision/ParticleCreationFunc.H" #include "Utils/TextMsg.H" @@ -67,6 +69,12 @@ CollisionHandler::CollisionHandler(MultiParticleContainer const * const mypc) collision_names[i], mypc ); } + else if (type == "bremsstrahlung") { + allcollisions[i] = + std::make_unique>( + collision_names[i], mypc + ); + } else{ WARPX_ABORT_WITH_MESSAGE("Unknown collision type."); } diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H index 7e3c89228ea..1423194bd72 100644 --- a/Source/Utils/ParticleUtils.H +++ b/Source/Utils/ParticleUtils.H @@ -124,6 +124,41 @@ namespace ParticleUtils { - gamma_V * Vz * gamma_u; } + /** + * \brief Perform a Lorentz transformation of the given velocity + * to a frame moving with gamma*velocity (Ux, Uy, Uz) relative to the present one. + * + * @param[in,out] ux,uy,uz components of velocity vector in the current + frame - importantly these quantities are gamma * velocity + * @param[in] Ux,Uy,Uz velocity of the new frame relative to the current one, + importantly these quantities are gamma * velocity + */ + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void doLorentzTransformWithU ( amrex::ParticleReal& ux, amrex::ParticleReal& uy, + amrex::ParticleReal& uz, + amrex::ParticleReal const Ux, amrex::ParticleReal const Uy, + amrex::ParticleReal const Uz ) + { + using namespace amrex::literals; + + constexpr auto c2 = PhysConst::c * PhysConst::c; + + amrex::ParticleReal const Usq = Ux*Ux + Uy*Uy + Uz*Uz; + amrex::ParticleReal const U = std::sqrt(Usq); + if (U > 0._prt) { + amrex::ParticleReal const gamma = std::sqrt(1._prt + Usq/c2); + amrex::ParticleReal const Ux_hat = Ux/U; + amrex::ParticleReal const Uy_hat = Uy/U; + amrex::ParticleReal const Uz_hat = Uz/U; + amrex::ParticleReal const P0 = std::sqrt(1._prt + (ux*ux + uy*uy + uz*uz)/c2); + amrex::ParticleReal const udotv = ux*Ux_hat + uy*Uy_hat + uz*Uz_hat; + ux = ux + (gamma - 1._prt)*udotv*Ux_hat - Ux*P0; + uy = uy + (gamma - 1._prt)*udotv*Uy_hat - Uy*P0; + uz = uz + (gamma - 1._prt)*udotv*Uz_hat - Uz*P0; + } + + } + /** * \brief Generate random unit vector in 3 dimensions * https://mathworld.wolfram.com/SpherePointPicking.html From 7f0f53ee0453242c0e3b82398830f4cf9c7fd591 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 3 Jan 2025 23:07:45 +0000 Subject: [PATCH 02/38] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 2 -- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp | 1 - 2 files changed, 3 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index a21ef31e72d..daef96bada6 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -473,5 +473,3 @@ private: }; #endif // WARPX_BREMSSTRAHLUNG_FUNC_H_ - - diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 5945127d013..2fc737316ba 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -63,4 +63,3 @@ BremsstrahlungFunc::UploadCrossSection (int Z) } } - From 03d96e05272936ced926d5c8e85454836db7fe75 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 3 Jan 2025 15:23:53 -0800 Subject: [PATCH 03/38] Fix typo --- .../Collision/BinaryCollision/Bremsstrahlung/Make.package | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package index 7bc4c15591c..b0792b413c1 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/Make.package @@ -1,4 +1,4 @@ CEXE_sources += BremsstrahlungFunc.cpp CEXE_sources += PhotonCreationFunc.cpp -VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision/Bremmstrahlung +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision/BinaryCollision/Bremsstrahlung From edf04f90d4d01f026a703a2463bdffb66d816793 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 3 Jan 2025 15:26:21 -0800 Subject: [PATCH 04/38] Small cleanup --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index daef96bada6..5749be1c505 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -317,7 +317,7 @@ public: fmulti = fmulti/arg; if (fmulti < 1.0_prt) { fmulti = 1.0_prt; } } - //q12 = 1.0 - exp(-arg); + //q12 = 1.0 - exp(-arg) amrex::ParticleReal q12 = arg; // Get a random number From e1105c363f0e098b7bd509f604c2072743fe5e38 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 3 Jan 2025 16:00:43 -0800 Subject: [PATCH 05/38] Small fixes --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 5749be1c505..ba85cadbf1c 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -203,7 +203,7 @@ public: } // find lo-index for koT1 cutoff (will typically be 1) - koT1_cut = std::max(KEcut_eV/KErel_eV, 1.0e-4); + koT1_cut = std::max(KEcut_eV/KErel_eV, 1.0e-4_prt); i0_cut = 0; for (int i=1; i < nkoT1; i++) { if (m_koT1_grid[i] > koT1_cut) { @@ -309,10 +309,12 @@ public: amrex::ParticleReal const beta1 = std::sqrt(gbp1sq_rel)/gammap1_rel; // magnitude electron beta amrex::ParticleReal arg = fmulti*beta1*PhysConst::c*sigma_total*n2*dt*gamma_factor; if (arg > 1.0_prt) { +#ifndef AMREX_USE_GPU amrex::Print() << "Notice: arg = " << arg << " in interSpeciesBremsstrahlung" << "\n"; amrex::Print() << " beta1 = " << beta1 << "\n"; amrex::Print() << " n2 = " << n2 << "\n"; amrex::Print() << " sigma_total = " << sigma_total << "\n"; +#endif arg = 1.0_prt; fmulti = fmulti/arg; if (fmulti < 1.0_prt) { fmulti = 1.0_prt; } From 2926baa2dfd0856f3c29f4f655e66d9f18c01c36 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 3 Jan 2025 16:27:43 -0800 Subject: [PATCH 06/38] More cleanup, including const's and prt's --- .../Bremsstrahlung/BremsstrahlungFunc.H | 60 +++++++++---------- 1 file changed, 30 insertions(+), 30 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index ba85cadbf1c..84f42e8787d 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -54,7 +54,7 @@ public: * @param[in] mypc pointer to the MultiParticleContainer * @param[in] isSameSpecies whether the two colliding species are the same (should always be false) */ - BremsstrahlungFunc (std::string const& collision_name, MultiParticleContainer const * const mypc, + BremsstrahlungFunc (std::string const& collision_name, MultiParticleContainer const * mypc, bool /*isSameSpecies*/); struct Executor { @@ -85,7 +85,7 @@ public: index_type const I2s, index_type const I2e, index_type const* AMREX_RESTRICT I1, index_type const* AMREX_RESTRICT I2, - const SoaData_type& soa_1, const SoaData_type& soa_2, + SoaData_type const& soa_1, const SoaData_type& soa_2, GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, amrex::ParticleReal const n1, amrex::ParticleReal const n2, amrex::ParticleReal const /*T1*/, amrex::ParticleReal const /*T2*/, @@ -109,12 +109,12 @@ public: amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; // Number of macroparticles of each species - const index_type NI1 = I1e - I1s; - const index_type NI2 = I2e - I2s; + index_type const NI1 = I1e - I1s; + index_type const NI2 = I2e - I2s; // Only loop over the number of particles in the first species (the electrons) - const index_type max_N = NI1; - const index_type min_N = NI2; + index_type const max_N = NI1; + index_type const min_N = NI2; index_type pair_index = cell_start_pair + coll_idx; @@ -124,7 +124,7 @@ public: #endif index_type i1 = I1s + coll_idx; - index_type i2 = I2s + coll_idx; + index_type const i2 = I2s + coll_idx; // we will start from collision number = coll_idx and then add // stride (smaller set size) until we do all collisions (larger set size) // If there are more atoms than electrons, this loop executes once @@ -161,12 +161,12 @@ public: u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta); #endif - } + p_pair_indices_1[pair_index] = I1[i1]; + p_pair_indices_2[pair_index] = I2[i2]; + i1 += min_N; + pair_index += min_N; - p_pair_indices_1[pair_index] = I1[i1]; - p_pair_indices_2[pair_index] = I2[i2]; - i1 += min_N; - pair_index += min_N; + } } static int constexpr nkoT1 = 17; @@ -221,20 +221,20 @@ public: int index; if (KErel_eV >= m_KEgrid_eV[nkoT1-1]) { index = nkoT1 - 2; - w0 = 0._prt; - w1 = 1._prt; + w0 = 0.0_prt; + w1 = 1.0_prt; } else { // find index for electron energy interpolation index = amrex::bisect(m_KEgrid_eV.data(), 0, nKE, KErel_eV); // compute interpolation weights for k*dsigma/dk w0 = (m_KEgrid_eV[index+1] - KErel_eV)/(m_KEgrid_eV[index+1] - m_KEgrid_eV[index]); - w1 = 1.0 - w0; + w1 = 1.0_prt - w0; } // kdsigdk at the cutoff amrex::ParticleReal const w00 = (m_koT1_grid[i0_cut+1] - koT1_cut)/(m_koT1_grid[i0_cut+1] - m_koT1_grid[i0_cut]); - amrex::ParticleReal const w01 = 1.0 - w00; + amrex::ParticleReal const w01 = 1.0_prt - w00; amrex::ParticleReal const kdsigdk_cut = ( w00*(w0*m_kdsigdk[index][i0_cut] + w1*m_kdsigdk[index+1][i0_cut]) + w01*(w0*m_kdsigdk[index][i0_cut+1] + w1*m_kdsigdk[index+1][i0_cut+1])); @@ -246,7 +246,7 @@ public: amrex::ParticleReal const this_k = (m_koT1_grid[i] + m_koT1_grid_im1)*0.5_prt; amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + w1*m_kdsigdk[index+1][i]; // using centered k here mitigates divergece effect for k->0 - amrex::ParticleReal const this_dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5/this_k; + amrex::ParticleReal const this_dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/this_k; sigmaC[i] = sigmaC[i-1] + this_dsigdk*this_dk; m_koT1_grid_im1 = m_koT1_grid[i]; kdsigdk_im1 = kdsigdk_i; @@ -279,10 +279,10 @@ public: amrex::ParticleReal const m_e_eV = m1*PhysConst::c*PhysConst::c/PhysConst::q_e; // 0.511e6 // compute gamma for particles 1 (electron) and 2 (ion) - amrex::ParticleReal gbp1sq = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; - amrex::ParticleReal gbp2sq = (u2x*u2x + u2y*u2y + u2z*u2z)/c2; - amrex::ParticleReal gammap1 = std::sqrt(1.0_prt + gbp1sq); - amrex::ParticleReal gammap2 = std::sqrt(1.0_prt + gbp2sq); + amrex::ParticleReal const gbp1sq = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; + amrex::ParticleReal const gbp2sq = (u2x*u2x + u2y*u2y + u2z*u2z)/c2; + amrex::ParticleReal const gammap1 = std::sqrt(1.0_prt + gbp1sq); + amrex::ParticleReal const gammap2 = std::sqrt(1.0_prt + gbp2sq); // transform electron to frame where ion is at rest amrex::ParticleReal u1x_rel = u1x; @@ -291,9 +291,9 @@ public: ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, u2x, u2y, u2z); // compute electron KE in frame where ion is at rest - amrex::ParticleReal gbp1sq_rel = (u1x_rel*u1x_rel + u1y_rel*u1y_rel + u1z_rel*u1z_rel)/c2; - amrex::ParticleReal gammap1_rel = std::sqrt(1.0_prt + gbp1sq_rel); - amrex::ParticleReal KE_eV = m_e_eV*gbp1sq_rel/(1.0_prt + gammap1_rel); + amrex::ParticleReal const gbp1sq_rel = (u1x_rel*u1x_rel + u1y_rel*u1y_rel + u1z_rel*u1z_rel)/c2; + amrex::ParticleReal const gammap1_rel = std::sqrt(1.0_prt + gbp1sq_rel); + amrex::ParticleReal const KE_eV = m_e_eV*gbp1sq_rel/(1.0_prt + gammap1_rel); amrex::GpuArray sigmaC; sigmaC.fill(0._prt); @@ -320,7 +320,7 @@ public: if (fmulti < 1.0_prt) { fmulti = 1.0_prt; } } //q12 = 1.0 - exp(-arg) - amrex::ParticleReal q12 = arg; + amrex::ParticleReal const q12 = arg; // Get a random number amrex::ParticleReal const random_number = amrex::Random(engine); @@ -331,14 +331,14 @@ public: amrex::ParticleReal const w_photon = w1/fmulti; // compute lab-frame electron kinetic energy before emission - amrex::ParticleReal KE_before = m_e_eV*gbp1sq/(1.0_prt + gammap1); + amrex::ParticleReal const KE_before = m_e_eV*gbp1sq/(1.0_prt + gammap1); // get energy of photon (hbar*omega) amrex::ParticleReal const random_number2 = amrex::Random(engine); int const index = amrex::bisect(sigmaC.data(), i0_cut, nkoT1, random_number2); amrex::ParticleReal const w0 = (random_number2 - sigmaC[index])/(sigmaC[index+1] - sigmaC[index]); amrex::ParticleReal const Epi = (index == i0_cut ? koT1_cut : m_koT1_grid[index]); - amrex::ParticleReal const EphotonoT1 = (1._prt - w0)*Epi + w0*m_koT1_grid[index+1]; + amrex::ParticleReal const EphotonoT1 = (1.0_prt - w0)*Epi + w0*m_koT1_grid[index+1]; amrex::ParticleReal const Ephoton_eV = EphotonoT1*KE_eV; // adjust electron gamma and beta @@ -358,9 +358,9 @@ public: u1z = u1z_rel; // compute electron kinetic energy after emission - gbp1sq = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; - gammap1 = std::sqrt(1.0_prt + gbp1sq); - amrex::ParticleReal const KE_after = m_e_eV*gbp1sq/(1.0_prt + gammap1); + amrex::ParticleReal const gbp1sq_prime = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; + amrex::ParticleReal const gammap1_prime = std::sqrt(1.0_prt + gbp1sq_prime); + amrex::ParticleReal const KE_after = m_e_eV*gbp1sq_prime/(1.0_prt + gammap1_prime); // update energy lost to bremsstrahlung amrex::ParticleReal const deltaE_eV = KE_after - KE_before; From 0ab5247b4a656c382ce127cf6ef3228169890df6 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 6 Jan 2025 09:01:01 -0800 Subject: [PATCH 07/38] Small fixes --- .../Bremsstrahlung/BremsstrahlungFunc.H | 5 ++--- .../Bremsstrahlung/BremsstrahlungFunc.cpp | 15 +++++++++------ 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 84f42e8787d..d3f7d0208a9 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -184,7 +184,7 @@ public: */ AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal - CalculateCrossSection (amrex::ParticleReal KErel, + CalculateCrossSection (amrex::ParticleReal KErel_eV, amrex::ParticleReal n1, amrex::ParticleReal m1, amrex::GpuArray & sigmaC, @@ -193,7 +193,6 @@ public: { using namespace amrex::literals; - amrex::ParticleReal const KErel_eV = KErel/PhysConst::q_e; amrex::ParticleReal const wpe = PhysConst::q_e*std::sqrt(n1/(m1*PhysConst::ep0)); amrex::ParticleReal const KEcut_eV = PhysConst::hbar*wpe/PhysConst::q_e; @@ -219,7 +218,7 @@ public: // kdsigdk is linearly weighted in electron energy amrex::ParticleReal w0, w1; int index; - if (KErel_eV >= m_KEgrid_eV[nkoT1-1]) { + if (KErel_eV >= m_KEgrid_eV[nKE-1]) { index = nkoT1 - 2; w0 = 0.0_prt; w1 = 1.0_prt; diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 2fc737316ba..d5c02826a3c 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -13,8 +13,11 @@ using namespace amrex::literals; BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, MultiParticleContainer const * const mypc, - bool /*isSameSpecies*/) + bool isSameSpecies) { + WARPX_ALWAYS_ASSERT_WITH_MESSAGE(!isSameSpecies, + "BremsstrahlungFunc: The two colliding species must be different"); + const amrex::ParmParse pp_collision_name(collision_name); // Read in the number of electrons on the target @@ -26,12 +29,12 @@ BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, Multi auto& product_species = mypc->GetParticleContainerFromName(product_species_name); WARPX_ALWAYS_ASSERT_WITH_MESSAGE(product_species.AmIA(), - "The product species must be photons"); + "BremsstrahlungFunc: The product species must be photons"); amrex::ParticleReal multiplier = 1._prt; - pp_collision_name.get("multiplier", multiplier); + pp_collision_name.query("multiplier", multiplier); WARPX_ALWAYS_ASSERT_WITH_MESSAGE(multiplier >= 1., - "The multiplier must be greater than or equal to one"); + "BremsstrahlungFunc: The multiplier must be greater than or equal to one"); m_exe.m_multiplier = multiplier; // Fill in the m_kdsigdk array @@ -51,13 +54,13 @@ BremsstrahlungFunc::UploadCrossSection (int Z) std::vector> & kdsigdk = m_kdsigdk_map.at(2); // Convert Seltzer and Berger energy-weighted differential cross section to units of [m^2] - for (int iee=0; iee < m_exe.nKE; iee++) { + for (int iee=0; iee < Executor::nKE; iee++) { amrex::ParticleReal const E = m_exe.m_KEgrid_eV[iee]/m_e_eV; amrex::ParticleReal const gamma = 1.0_rt + E; /* betaSq = 1.0 - 1.0/gamma/gamma */ amrex::ParticleReal const betaSq = (E*E + 2._rt*E)/gamma/gamma; amrex::ParticleReal const scale_factor = 1.0e-31_rt*Z*Z/betaSq; - for (int iep=0; iep < m_exe.nkoT1; iep++) { + for (int iep=0; iep < Executor::nkoT1; iep++) { m_exe.m_kdsigdk[iee][iep] = kdsigdk[iee][iep]*scale_factor; } } From 8a147cb145cd47eaf4e88448c62d9d6bc6aaec52 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 6 Jan 2025 13:58:11 -0800 Subject: [PATCH 08/38] Bug fix --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index d5c02826a3c..1d7d34b8204 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -51,7 +51,7 @@ BremsstrahlungFunc::UploadCrossSection (int Z) constexpr auto m_e_eV = PhysConst::m_e*PhysConst::c*PhysConst::c/PhysConst::q_e; // 0.511e6 - std::vector> & kdsigdk = m_kdsigdk_map.at(2); + std::vector> & kdsigdk = m_kdsigdk_map.at(Z); // Convert Seltzer and Berger energy-weighted differential cross section to units of [m^2] for (int iee=0; iee < Executor::nKE; iee++) { From d3b22d7da9e7458e616ce66219df4f79e0e97688 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 6 Jan 2025 16:34:56 -0800 Subject: [PATCH 09/38] Fix a variable name --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index d3f7d0208a9..ccd9e2066ac 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -238,16 +238,16 @@ public: + w01*(w0*m_kdsigdk[index][i0_cut+1] + w1*m_kdsigdk[index+1][i0_cut+1])); // create cumulative distribution using trapezoidal rule - amrex::ParticleReal m_koT1_grid_im1 = koT1_cut; + amrex::ParticleReal koT1_grid_im1 = koT1_cut; amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; for (int i=i0_cut+1; i < nkoT1; i++) { - amrex::ParticleReal const this_dk = (m_koT1_grid[i] - m_koT1_grid_im1); - amrex::ParticleReal const this_k = (m_koT1_grid[i] + m_koT1_grid_im1)*0.5_prt; + amrex::ParticleReal const this_dk = (m_koT1_grid[i] - koT1_grid_im1); + amrex::ParticleReal const this_k = (m_koT1_grid[i] + koT1_grid_im1)*0.5_prt; amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + w1*m_kdsigdk[index+1][i]; // using centered k here mitigates divergece effect for k->0 amrex::ParticleReal const this_dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/this_k; sigmaC[i] = sigmaC[i-1] + this_dsigdk*this_dk; - m_koT1_grid_im1 = m_koT1_grid[i]; + koT1_grid_im1 = m_koT1_grid[i]; kdsigdk_im1 = kdsigdk_i; } amrex::ParticleReal const sigma_total = sigmaC[nkoT1-1]; // total cross section [m^2] From 2dcd95adde0f1ba673d9db61be77791cb36b217e Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 6 Jan 2025 16:49:29 -0800 Subject: [PATCH 10/38] Initialize p_mask to zero Since the number of collision is the number of electrons there are some pairs the are not done when the number of species 2 is greater than the number of electrons. --- Source/Particles/Collision/BinaryCollision/BinaryCollision.H | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H index b4eaadeee76..dad819bd12e 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollision.H +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollision.H @@ -581,7 +581,7 @@ public: index_type* AMREX_RESTRICT p_coll_offsets = coll_offsets.dataPtr(); // mask: equal to 1 if particle creation occurs for a given pair, 0 otherwise - amrex::Gpu::DeviceVector mask(n_total_pairs); + amrex::Gpu::DeviceVector mask(n_total_pairs, index_type(0)); index_type* AMREX_RESTRICT p_mask = mask.dataPtr(); // Will be filled with the index of the first particle of a given pair amrex::Gpu::DeviceVector pair_indices_1(n_total_pairs); From c50b855c732855b9c1ea77775f30b66bf355d64e Mon Sep 17 00:00:00 2001 From: David Grote Date: Wed, 8 Jan 2025 14:58:57 -0800 Subject: [PATCH 11/38] Fix the energy scaling of the generated photons --- .../BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H index ac9c213b104..77e2bb42b1e 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H @@ -168,10 +168,11 @@ public: auto& uz1 = soa_products_data[0].m_rdata[PIdx::uz][product_index]; // Normalize out the electron velocity and multiply by the photon momentum, E/c + // Also, the photon momentum is normalized by m_e auto u1 = std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1); - ux1 *= (p_product_data[i]/PhysConst::c)/u1; - uy1 *= (p_product_data[i]/PhysConst::c)/u1; - uz1 *= (p_product_data[i]/PhysConst::c)/u1; + ux1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; + uy1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; + uz1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; } }); From 1bb5422975c9e9d26e821d1495261cef82ded68e Mon Sep 17 00:00:00 2001 From: David Grote Date: Wed, 8 Jan 2025 15:00:47 -0800 Subject: [PATCH 12/38] Improvements to the algorithm Removed the use of the temporary sigmaC array. Now use quadratic weighting for the photon energy. --- .../Bremsstrahlung/BremsstrahlungFunc.H | 134 +++++++++++++----- 1 file changed, 96 insertions(+), 38 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index ccd9e2066ac..1a3e24f3b53 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -98,12 +98,12 @@ public: amrex::ParticleReal* AMREX_RESTRICT p_product_data, amrex::RandomEngine const& engine) const { - amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w]; + amrex::ParticleReal * const AMREX_RESTRICT weight1 = soa_1.m_rdata[PIdx::w]; amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux]; amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz]; - amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w]; + amrex::ParticleReal * const AMREX_RESTRICT weight2 = soa_2.m_rdata[PIdx::w]; amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux]; amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy]; amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz]; @@ -150,7 +150,7 @@ public: BremsstrahlungEvent( u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], - m1, w1[ I1[i1] ], w2[ I2[i2] ], + m1, weight1[ I1[i1] ], weight2[ I2[i2] ], n1, n2, dt, static_cast(pair_index), p_mask, p_pair_reaction_weight, p_product_data, engine); @@ -175,11 +175,15 @@ public: /* \brief Calculate the cross section given the electron energy. * The differential cross section is cut off below the plasma frequency since the * low energy photons would all be immediately absorbed by the surrounding plasma. - * The normalized sigmaC is returned to allow calculating the photon energy. * * \param[in] E electron energy in Joules * \param[in] n1 the electron number density - * \param[out] sigmaC the normalized differential cross section + * \param[in] m1 the mass of the electrons + * \param[in] index the index for the electron energy grid + * \param[in] i0_cut the index below the cut off for the photon energy grid + * \param[in] koT1_cut the k of the photon energy cut off + * \param[in] kdsigdk_cut the ksigdk at the cut off + * \param[in] w0 the grid weighting at the cut off * \result the cross section */ AMREX_GPU_HOST_DEVICE AMREX_INLINE @@ -187,9 +191,10 @@ public: CalculateCrossSection (amrex::ParticleReal KErel_eV, amrex::ParticleReal n1, amrex::ParticleReal m1, - amrex::GpuArray & sigmaC, + int & index, int & i0_cut, - amrex::ParticleReal & koT1_cut) const + amrex::ParticleReal & koT1_cut, amrex::ParticleReal & kdsigdk_cut, + amrex::ParticleReal & w0) const { using namespace amrex::literals; @@ -216,55 +221,110 @@ public: } // kdsigdk is linearly weighted in electron energy - amrex::ParticleReal w0, w1; - int index; if (KErel_eV >= m_KEgrid_eV[nKE-1]) { index = nkoT1 - 2; w0 = 0.0_prt; - w1 = 1.0_prt; } else { // find index for electron energy interpolation index = amrex::bisect(m_KEgrid_eV.data(), 0, nKE, KErel_eV); // compute interpolation weights for k*dsigma/dk w0 = (m_KEgrid_eV[index+1] - KErel_eV)/(m_KEgrid_eV[index+1] - m_KEgrid_eV[index]); - w1 = 1.0_prt - w0; } // kdsigdk at the cutoff amrex::ParticleReal const w00 = (m_koT1_grid[i0_cut+1] - koT1_cut)/(m_koT1_grid[i0_cut+1] - m_koT1_grid[i0_cut]); amrex::ParticleReal const w01 = 1.0_prt - w00; - amrex::ParticleReal const kdsigdk_cut = ( w00*(w0*m_kdsigdk[index][i0_cut] + w1*m_kdsigdk[index+1][i0_cut]) - + w01*(w0*m_kdsigdk[index][i0_cut+1] + w1*m_kdsigdk[index+1][i0_cut+1])); + kdsigdk_cut = ( w00*(w0*m_kdsigdk[index][i0_cut ] + (1.0_prt - w0)*m_kdsigdk[index+1][i0_cut]) + + w01*(w0*m_kdsigdk[index][i0_cut+1] + (1.0_prt - w0)*m_kdsigdk[index+1][i0_cut+1])); + // Scale kdsigdk_cut by the centered k to mitigate divergence effect for k->0 + kdsigdk_cut *= koT1_cut/((m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt); - // create cumulative distribution using trapezoidal rule + // Integrate the distribution using trapezoidal rule amrex::ParticleReal koT1_grid_im1 = koT1_cut; amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; + amrex::ParticleReal sigma = 0.0_prt; for (int i=i0_cut+1; i < nkoT1; i++) { - amrex::ParticleReal const this_dk = (m_koT1_grid[i] - koT1_grid_im1); - amrex::ParticleReal const this_k = (m_koT1_grid[i] + koT1_grid_im1)*0.5_prt; - amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + w1*m_kdsigdk[index+1][i]; - // using centered k here mitigates divergece effect for k->0 - amrex::ParticleReal const this_dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/this_k; - sigmaC[i] = sigmaC[i-1] + this_dsigdk*this_dk; + amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1); + amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i]; + amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt; + sigma = sigma + dsigdk*dk; koT1_grid_im1 = m_koT1_grid[i]; kdsigdk_im1 = kdsigdk_i; } - amrex::ParticleReal const sigma_total = sigmaC[nkoT1-1]; // total cross section [m^2] + amrex::ParticleReal const sigma_total = sigma; // total cross section [m^2] - for (int i=i0_cut+1; i < nkoT1; i++) { - sigmaC[i] /= sigma_total; + return sigma_total; + } + + /* \brief Generate the photon energy from the differential cross section + * + * \param[in] index the index for the electron energy grid + * \param[in] i0_cut the index below the cut off for the photon energy grid + * \param[in] koT1_cut the k of the photon energy cut off + * \param[in] kdsigdk_cut the ksigdk at the cut off + * \param[in] w0 the grid weighting at the cut off + * \param[in] sigma_total the total cross section + * \param[in] random_number the uniformly spaced random number + * \result the photon energy + */ + AMREX_GPU_HOST_DEVICE AMREX_INLINE + amrex::ParticleReal + Photon_energy(int index, int i0_cut, amrex::ParticleReal koT1_cut, amrex::ParticleReal kdsigdk_cut, + amrex::ParticleReal w0, amrex::ParticleReal sigma_total, amrex::ParticleReal random_number) const + { + using namespace amrex::literals; + amrex::ParticleReal koT1_grid_im1 = koT1_cut; + amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; + amrex::ParticleReal sigma = 0._prt; + amrex::ParticleReal kdsigdk_i; + amrex::ParticleReal dk; + int i; + for (i=i0_cut+1; i < nkoT1; i++) { + dk = (m_koT1_grid[i] - koT1_grid_im1); + kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i]; + amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt; + amrex::ParticleReal const next_sigma = sigma + dsigdk*dk/sigma_total; + if (random_number > next_sigma) { + sigma = next_sigma; + koT1_grid_im1 = m_koT1_grid[i]; + kdsigdk_im1 = kdsigdk_i; + } else { + break; + } } - return sigma_total; + amrex::ParticleReal const nnorm_i = kdsigdk_im1/koT1_grid_im1/sigma_total; + amrex::ParticleReal const nnorm_ip1 = kdsigdk_i/m_koT1_grid[i]/sigma_total; + amrex::ParticleReal const result = ((std::sqrt(nnorm_i*nnorm_i + + 2._prt*(nnorm_ip1 - nnorm_i)*(random_number - sigma)/dk) + - nnorm_i)/(nnorm_ip1 - nnorm_i))*dk + koT1_grid_im1; + + return result; } + /* \brief Do the Bremsstrahlung calculation + * + * \param[in] u1x, u1y, u1z the gamma*velocity of the first species, the electrons + * \param[in] u2x, u2y, u2z the gamma*velocity of the second species, the atoms + * \param[in] m1 the mass of the electrons + * \param[in] weight1 the particle weight of the electrons + * \param[in] weight2 the particle weight of the target (unused) + * \param[in] n1 the number density of the electrons + * \param[in] n2 the number density of the target + * \param[in] dt the time step size + * \param[in] pair_index the index number of the pair colliding + * \param[out] p_mask flags whether the pair collided + * \param[out] p_pair_reaction_weight the particle weight of the generated photon + * \param[out] p_product_data the energy of the generated photon + * \param[in] engine the random number generating engine + */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void BremsstrahlungEvent (amrex::ParticleReal & u1x, amrex::ParticleReal & u1y, amrex::ParticleReal & u1z, amrex::ParticleReal const& u2x, amrex::ParticleReal const& u2y, amrex::ParticleReal const& u2z, amrex::ParticleReal m1, - amrex::ParticleReal w1, amrex::ParticleReal /*w2*/, + amrex::ParticleReal weight1, amrex::ParticleReal /*weight2*/, amrex::ParticleReal n1, amrex::ParticleReal n2, amrex::Real dt, int pair_index, index_type* AMREX_RESTRICT p_mask, @@ -294,11 +354,11 @@ public: amrex::ParticleReal const gammap1_rel = std::sqrt(1.0_prt + gbp1sq_rel); amrex::ParticleReal const KE_eV = m_e_eV*gbp1sq_rel/(1.0_prt + gammap1_rel); - amrex::GpuArray sigmaC; - sigmaC.fill(0._prt); - int i0_cut; - amrex::ParticleReal koT1_cut; - amrex::ParticleReal const sigma_total = CalculateCrossSection(KE_eV, n1, m1, sigmaC, i0_cut, koT1_cut); + int i0_cut, index; + amrex::ParticleReal koT1_cut, kdsigdk_cut; + amrex::ParticleReal w0; + amrex::ParticleReal const sigma_total = CalculateCrossSection(KE_eV, n1, m1, + index, i0_cut, koT1_cut, kdsigdk_cut, w0); if (sigma_total == 0._prt) { return; } @@ -327,21 +387,19 @@ public: if (random_number < q12) { // compute weight for photon - amrex::ParticleReal const w_photon = w1/fmulti; + amrex::ParticleReal const w_photon = weight1/fmulti; // compute lab-frame electron kinetic energy before emission amrex::ParticleReal const KE_before = m_e_eV*gbp1sq/(1.0_prt + gammap1); // get energy of photon (hbar*omega) amrex::ParticleReal const random_number2 = amrex::Random(engine); - int const index = amrex::bisect(sigmaC.data(), i0_cut, nkoT1, random_number2); - amrex::ParticleReal const w0 = (random_number2 - sigmaC[index])/(sigmaC[index+1] - sigmaC[index]); - amrex::ParticleReal const Epi = (index == i0_cut ? koT1_cut : m_koT1_grid[index]); - amrex::ParticleReal const EphotonoT1 = (1.0_prt - w0)*Epi + w0*m_koT1_grid[index+1]; + amrex::ParticleReal const EphotonoT1 = Photon_energy(index, i0_cut, koT1_cut, kdsigdk_cut, + w0, sigma_total, random_number2); amrex::ParticleReal const Ephoton_eV = EphotonoT1*KE_eV; // adjust electron gamma and beta - amrex::ParticleReal const KE_prime_eV = KE_eV - Ephoton_eV*w_photon/w1; + amrex::ParticleReal const KE_prime_eV = KE_eV - Ephoton_eV*w_photon/weight1; amrex::ParticleReal const gamma_prime = 1.0_prt + KE_prime_eV/m_e_eV; // rescale electron beta in frame where ion is at rest @@ -363,10 +421,10 @@ public: // update energy lost to bremsstrahlung amrex::ParticleReal const deltaE_eV = KE_after - KE_before; - /* m_deltaE_bremsstrahlung -= deltaE_eV*w1*PhysConst::q_e; // Joules */ + /* m_deltaE_bremsstrahlung -= deltaE_eV*weight1*PhysConst::q_e; // Joules */ // compute photon energy in lab frame - amrex::ParticleReal const Ephoton_lab = deltaE_eV*w1/w_photon*PhysConst::q_e; // Joules + amrex::ParticleReal const Ephoton_lab = deltaE_eV*weight1/w_photon*PhysConst::q_e; // Joules p_pair_reaction_weight[pair_index] = w_photon; p_product_data[pair_index] = Ephoton_lab; From 27fe428b64e25de180770abe9d977a01b71aafa3 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 13 Jan 2025 10:44:07 -0800 Subject: [PATCH 13/38] Add CI test --- Examples/Tests/collision/CMakeLists.txt | 10 ++ .../analysis_collision_1d_Bremsstrahlung.py | 86 ++++++++++++++ .../inputs_test_1d_collision_z_Bremsstrahlung | 112 ++++++++++++++++++ .../test_1d_collision_z_Bremsstrahlung.json | 34 ++++++ 4 files changed, 242 insertions(+) create mode 100755 Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py create mode 100644 Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung create mode 100644 Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json diff --git a/Examples/Tests/collision/CMakeLists.txt b/Examples/Tests/collision/CMakeLists.txt index 522dafbfbfb..9b6bff17d53 100644 --- a/Examples/Tests/collision/CMakeLists.txt +++ b/Examples/Tests/collision/CMakeLists.txt @@ -60,3 +60,13 @@ add_warpx_test( "analysis_default_regression.py --path diags/diag1000150 --skip-particles" # checksum OFF # dependency ) + +add_warpx_test( + test_1d_collision_z_Bremsstrahlung # name + 1 # dims + 1 # nprocs + inputs_test_1d_collision_z_Bremsstrahlung # inputs + "analysis_collision_1d_Bremsstrahlung.py diags/diag1000100" # analysis + "analysis_default_regression.py --path diags/diag1000100" # checksum + OFF # dependency +) diff --git a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py new file mode 100755 index 00000000000..69835cea60f --- /dev/null +++ b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py @@ -0,0 +1,86 @@ +#!/usr/bin/env python3 + +# Copyright 2025 David Grote +# +# +# This file is part of WarpX. +# +# License: BSD-3-Clause-LBNL +# +# This is a script that analyses the simulation results from the script `inputs_test_1d_collision_z_Bremsstrahlung`. +# run locally: python analysis_collision_1d_Bremsstrahlung.py diags/diag1000600/ +# +# This is a 1D Bremsstrahlung collision of electrons on Borons ions +# producing photons. This checks that the appropriate number of +# photons are created with the correct distribution. +# +import sys +import os + +import numpy as np +from scipy import constants + +# this will be the name of the plot file +last_fn = sys.argv[1] + +particle_energy = np.loadtxt(os.path.join("diags", "reducedfiles", "particle_energy.txt"), skiprows=1) +particle_number = np.loadtxt(os.path.join("diags", "reducedfiles", "particle_number.txt"), skiprows=1) +photon_energy = np.loadtxt(os.path.join("diags", "reducedfiles", "photon_energy.txt"), skiprows=1) + +total_energy = particle_energy[:,2] +electron_energy = particle_energy[:,3] +ion_energy = particle_energy[:,4] +photon_energy = particle_energy[:,5] + +energy_tolerance = 1.e-12 + +print(f"initial total energy = {total_energy[0]}") +print(f"final total energy = {total_energy[-1]}") + +dE_total = np.abs(total_energy[-1] - total_energy[0])/total_energy[0] +print(f"change in total energy = {dE_total}") +assert dE_total < energy_tolerance + +print(f"initial electron energy = {electron_energy[0]}") +print(f"final electron energy = {electron_energy[-1]}") +print(f"final photon energy = {photon_energy[-1]}") + +dE_electron_photon = np.abs(electron_energy[-1] + photon_energy[-1] - electron_energy[0])/electron_energy[0] +print(f"electron, photon energy change = {dE_electron_photon}") +assert dE_electron_photon < energy_tolerance + + +print() + +dt = 1.e-2*1.e-15 +Z = 5 +n_i = 5.47e31 +n_e = 5.47e30 +L = 1.e-6 # 1 micron +T1 = 1.e6 # 1 MeV +N_e = n_e*L # number of electrons +m_e_eV = constants.m_e*constants.c**2/constants.e +gamma = T1/m_e_eV + 1. +gamma_beta = np.sqrt(gamma**2 - 1.) +beta = gamma_beta/gamma + +phirad = 6.761 # from Seltzer and Berger for 1 MeV electron and Boron + + +dEdx_simulation = (particle_energy[0,3] - particle_energy[1:,3])/particle_energy[1:,0]/(beta*constants.c*dt)/constants.e/N_e + +Boron_weight = 20065.0*constants.m_e +r_e = 1./(4.*constants.pi*constants.epsilon_0)*(constants.e**2/(constants.m_e*constants.c**2)) +dEdx = n_i*constants.alpha*r_e**2*Z**2*(T1 + m_e_eV)*phirad + +print(f"dE/dx analytic = {dEdx}") +print(F"dE/dx simulated = {dEdx_simulation[-1]}") +print(f"dE/dx simulated/analytic = {dEdx_simulation[-1]/dEdx}") +assert np.abs(dEdx_simulation[-1]/dEdx - 1.) < 0.03 + +sigam_total = 1.698e-28 # Calculated from table with k cutoff=1.e-4 +N_photon = n_e*n_i*L*beta*constants.c*sigam_total*dt +print(f"New photons per step simulated/analytic = {particle_number[-1,-1]/(particle_energy[-1,0]*N_photon)}") +assert np.abs(particle_number[-1,-1]/(particle_energy[-1,0]*N_photon) - 1.) < 0.03 + + diff --git a/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung new file mode 100644 index 00000000000..ba2955a5f6d --- /dev/null +++ b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung @@ -0,0 +1,112 @@ +################################# +########## CONSTANTS ############ +################################# + +my_constants.n_e = 5.47e30 # electron density, m^-3 +my_constants.Np_e = 1024 # electrons per cell +my_constants.T1 = 1.e6 # 1 MeV, electron energy +my_constants.m_e_eV = m_e*clight**2/q_e # electron rest mass in eV +my_constants.gamma = T1/m_e_eV + 1. # electron gamma factor +my_constants.gamma_beta = sqrt(gamma**2 - 1.) + +my_constants.n_i = 5.47e31 # ion density, m^-3, for rho = 1000 g/cm^3 +my_constants.Np_i = 1024 # ions per cell +my_constants.T_i = 2.0 # ion temperature, eV + +my_constants.m_B11 = 20065.0*m_e # 11.01*amu/me - 5, Boron +my_constants.q_B11 = 5.*q_e + +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 100 +amr.n_cell = 100 +amr.max_level = 0 +amr.blocking_factor = 4 +geometry.dims = 1 +geometry.prob_lo = 0. +geometry.prob_hi = 1.e-6 + +################################# +###### Boundary Condition ####### +################################# +boundary.field_lo = periodic +boundary.field_hi = periodic + +################################# +############ NUMERICS ########### +################################# +warpx.serialize_initial_conditions = 1 +warpx.verbose = 1 +warpx.const_dt = 1.e-2*1.e-15 +warpx.use_filter = 0 + +# Do not evolve the E and B fields +algo.maxwell_solver = none + +# Order of particle shape factors +algo.particle_shape = 1 + +################################# +############ PLASMA ############# +################################# +particles.species_names = electrons ions photons +particles.photon_species = photons + +electrons.species_type = electron +electrons.do_not_deposit = 1 + +electrons.injection_style = "NUniformPerCell" +electrons.num_particles_per_cell_each_dim = Np_e +electrons.profile = "constant" +electrons.density = n_e +electrons.momentum_distribution_type = "constant" +electrons.uz = gamma_beta + +ions.charge = q_B11 +ions.mass = m_B11 +ions.do_not_deposit = 1 + +ions.injection_style = "NUniformPerCell" +ions.num_particles_per_cell_each_dim = Np_i +ions.profile = "constant" +ions.density = n_i +ions.momentum_distribution_type = "gaussian" +ions.ux_th = sqrt(T_i*q_e/m_B11)/clight +ions.uy_th = sqrt(T_i*q_e/m_B11)/clight +ions.uz_th = sqrt(T_i*q_e/m_B11)/clight + +photons.species_type = photon +photons.injection_style = none + +################################# +############ COLLISION ########## +################################# +collisions.collision_names = bremsstrahlung +bremsstrahlung.type = bremsstrahlung +bremsstrahlung.species = electrons ions +bremsstrahlung.Z = 5 +bremsstrahlung.product_species = photons +bremsstrahlung.multiplier = 100. + +# Diagnostics +diagnostics.diags_names = diag1 +diag1.intervals = 100 +diag1.diag_type = Full + +warpx.reduced_diags_names = particle_energy particle_number photon_energy +particle_energy.type = ParticleEnergy +particle_energy.intervals = 10 +particle_energy.precision = 18 +particle_number.type = ParticleNumber +particle_number.intervals = 10 +particle_number.precision = 18 +photon_energy.type = ParticleHistogram +photon_energy.species = photons +photon_energy.histogram_function(t,x,y,z,ux,uy,uz) = log10(sqrt(ux*ux+uy*uy+uz*uz)*clight*clight/q_e*m_e) +photon_energy.bin_number = 100 +photon_energy.bin_min = log10(10.) +photon_energy.bin_max = log10(T1) +photon_energy.intervals = 10 +photon_energy.precision = 18 + diff --git a/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json b/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json new file mode 100644 index 00000000000..d8edba3a5b4 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json @@ -0,0 +1,34 @@ +{ + "lev=0": { + "Bx": 0.0, + "By": 0.0, + "Bz": 0.0, + "Ex": 0.0, + "Ey": 0.0, + "Ez": 0.0, + "jx": 0.0, + "jy": 0.0, + "jz": 0.0 + }, + "electrons": { + "particle_momentum_x": 1.7785184549415273e-26, + "particle_momentum_y": 1.769780531165544e-26, + "particle_momentum_z": 7.780435420929519e-17, + "particle_position_x": 0.0511996620835328, + "particle_weight": 5.470000000000001e+24 + }, + "photons": { + "particle_momentum_x": 3.9236321845123015e-27, + "particle_momentum_y": 3.892475715179289e-27, + "particle_momentum_z": 1.2828655290730977e-18, + "particle_position_x": 0.013297153040958081, + "particle_weight": 1.4269541992187505e+22 + }, + "ions": { + "particle_momentum_x": 6.245889688274839e-18, + "particle_momentum_y": 6.2336526718397175e-18, + "particle_momentum_z": 6.258023561029987e-18, + "particle_position_x": 0.05120000226548876, + "particle_weight": 5.470000000000001e+25 + } +} From 3d262458f40265c4131cf5a122dcbddb4eb223aa Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 13 Jan 2025 10:58:29 -0800 Subject: [PATCH 14/38] Make creation of photons optional --- .../Bremsstrahlung/BremsstrahlungFunc.H | 12 ++++++++---- .../Bremsstrahlung/BremsstrahlungFunc.cpp | 4 ++++ 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 1a3e24f3b53..584763e4d80 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -423,11 +423,14 @@ public: amrex::ParticleReal const deltaE_eV = KE_after - KE_before; /* m_deltaE_bremsstrahlung -= deltaE_eV*weight1*PhysConst::q_e; // Joules */ - // compute photon energy in lab frame - amrex::ParticleReal const Ephoton_lab = deltaE_eV*weight1/w_photon*PhysConst::q_e; // Joules + if (m_create_photons) { + // compute photon energy in lab frame + amrex::ParticleReal const Ephoton_lab = deltaE_eV*weight1/w_photon*PhysConst::q_e; // Joules - p_pair_reaction_weight[pair_index] = w_photon; - p_product_data[pair_index] = Ephoton_lab; + p_mask[pair_index] = true; + p_pair_reaction_weight[pair_index] = w_photon; + p_product_data[pair_index] = Ephoton_lab; + } } @@ -435,6 +438,7 @@ public: bool m_computeSpeciesDensities = true; bool m_computeSpeciesTemperatures = false; + bool m_create_photons = true; bool m_need_product_data = true; amrex::ParticleReal m_multiplier; diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 1d7d34b8204..58e3ef27048 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -31,6 +31,10 @@ BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, Multi WARPX_ALWAYS_ASSERT_WITH_MESSAGE(product_species.AmIA(), "BremsstrahlungFunc: The product species must be photons"); + bool create_photons = true; + pp_collision_name.query("create_photons", create_photons); + m_exe.m_create_photons = create_photons + amrex::ParticleReal multiplier = 1._prt; pp_collision_name.query("multiplier", multiplier); WARPX_ALWAYS_ASSERT_WITH_MESSAGE(multiplier >= 1., From 43939ca22c8dcf0c1381cc27b542de29090d7b1d Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 13 Jan 2025 11:33:56 -0800 Subject: [PATCH 15/38] Add documentation --- Docs/source/usage/parameters.rst | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 31f0e06ab5b..a91cb03c815 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -2059,17 +2059,21 @@ Details about the collision models can be found in the :ref:`theory section `__. * ``.species`` (`strings`) - If using ``dsmc``, ``pairwisecoulomb`` or ``nuclearfusion``, this should be the name(s) of the species, + If using ``dsmc``, ``pairwisecoulomb``, ``nuclearfusion``, or ``bremsstrahlung``, this should be the name(s) of the species, between which the collision will be considered. (Provide only one name for intra-species collisions.) + With ``bremsstrahlung``, the electron species must be given first, followed by the target species. If using ``background_mcc`` or ``background_stopping`` type this should be the name of the species for which collisions with a background will be included. In this case, only one species name should be given. * ``.product_species`` (`strings`) - Only for ``nuclearfusion``. The name(s) of the species in which to add + Only for ``nuclearfusion`` and ``bremsstrahlung``. The name(s) of the species in which to add the new macroparticles created by the reaction. + For ``bremsstrahlung``, the product species must be of type photon. * ``.ndt`` (`int`) optional Execute collision every # time steps. The default value is 1. @@ -2191,6 +2195,18 @@ Details about the collision models can be found in the :ref:`theory section .Z`` (`integer`) + Only for ``bremsstrahlung``. The atomic number of the target ion species. + Currently, only the values 1, 2, 5, 6 are supported. + +* ``.multiplier`` (`float`) + Only for ``bremsstrahlung``. Multiplier for the collision probability. + Any resulting photons will have the electron weight divided the multiplier. + The default is 1. This must be greater than or equal to 1. + +* ``.create_photons`` (`integer`) + Only for ``bremsstrahlung``. Whether photons will be created, defaults to 1 (true). + .. _running-cpp-parameters-numerics: Numerics and algorithms From fd4a0bece4f2006eff70942d2877a60afa75c688 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 13 Jan 2025 19:42:31 +0000 Subject: [PATCH 16/38] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../analysis_collision_1d_Bremsstrahlung.py | 81 ++++++++++++------- .../inputs_test_1d_collision_z_Bremsstrahlung | 1 - 2 files changed, 51 insertions(+), 31 deletions(-) diff --git a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py index 69835cea60f..d3242d6c91e 100755 --- a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py +++ b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py @@ -14,8 +14,8 @@ # producing photons. This checks that the appropriate number of # photons are created with the correct distribution. # -import sys import os +import sys import numpy as np from scipy import constants @@ -23,21 +23,27 @@ # this will be the name of the plot file last_fn = sys.argv[1] -particle_energy = np.loadtxt(os.path.join("diags", "reducedfiles", "particle_energy.txt"), skiprows=1) -particle_number = np.loadtxt(os.path.join("diags", "reducedfiles", "particle_number.txt"), skiprows=1) -photon_energy = np.loadtxt(os.path.join("diags", "reducedfiles", "photon_energy.txt"), skiprows=1) +particle_energy = np.loadtxt( + os.path.join("diags", "reducedfiles", "particle_energy.txt"), skiprows=1 +) +particle_number = np.loadtxt( + os.path.join("diags", "reducedfiles", "particle_number.txt"), skiprows=1 +) +photon_energy = np.loadtxt( + os.path.join("diags", "reducedfiles", "photon_energy.txt"), skiprows=1 +) -total_energy = particle_energy[:,2] -electron_energy = particle_energy[:,3] -ion_energy = particle_energy[:,4] -photon_energy = particle_energy[:,5] +total_energy = particle_energy[:, 2] +electron_energy = particle_energy[:, 3] +ion_energy = particle_energy[:, 4] +photon_energy = particle_energy[:, 5] -energy_tolerance = 1.e-12 +energy_tolerance = 1.0e-12 print(f"initial total energy = {total_energy[0]}") print(f"final total energy = {total_energy[-1]}") -dE_total = np.abs(total_energy[-1] - total_energy[0])/total_energy[0] +dE_total = np.abs(total_energy[-1] - total_energy[0]) / total_energy[0] print(f"change in total energy = {dE_total}") assert dE_total < energy_tolerance @@ -45,42 +51,57 @@ print(f"final electron energy = {electron_energy[-1]}") print(f"final photon energy = {photon_energy[-1]}") -dE_electron_photon = np.abs(electron_energy[-1] + photon_energy[-1] - electron_energy[0])/electron_energy[0] +dE_electron_photon = ( + np.abs(electron_energy[-1] + photon_energy[-1] - electron_energy[0]) + / electron_energy[0] +) print(f"electron, photon energy change = {dE_electron_photon}") assert dE_electron_photon < energy_tolerance print() -dt = 1.e-2*1.e-15 +dt = 1.0e-2 * 1.0e-15 Z = 5 n_i = 5.47e31 n_e = 5.47e30 -L = 1.e-6 # 1 micron -T1 = 1.e6 # 1 MeV -N_e = n_e*L # number of electrons -m_e_eV = constants.m_e*constants.c**2/constants.e -gamma = T1/m_e_eV + 1. -gamma_beta = np.sqrt(gamma**2 - 1.) -beta = gamma_beta/gamma +L = 1.0e-6 # 1 micron +T1 = 1.0e6 # 1 MeV +N_e = n_e * L # number of electrons +m_e_eV = constants.m_e * constants.c**2 / constants.e +gamma = T1 / m_e_eV + 1.0 +gamma_beta = np.sqrt(gamma**2 - 1.0) +beta = gamma_beta / gamma phirad = 6.761 # from Seltzer and Berger for 1 MeV electron and Boron -dEdx_simulation = (particle_energy[0,3] - particle_energy[1:,3])/particle_energy[1:,0]/(beta*constants.c*dt)/constants.e/N_e +dEdx_simulation = ( + (particle_energy[0, 3] - particle_energy[1:, 3]) + / particle_energy[1:, 0] + / (beta * constants.c * dt) + / constants.e + / N_e +) -Boron_weight = 20065.0*constants.m_e -r_e = 1./(4.*constants.pi*constants.epsilon_0)*(constants.e**2/(constants.m_e*constants.c**2)) -dEdx = n_i*constants.alpha*r_e**2*Z**2*(T1 + m_e_eV)*phirad +Boron_weight = 20065.0 * constants.m_e +r_e = ( + 1.0 + / (4.0 * constants.pi * constants.epsilon_0) + * (constants.e**2 / (constants.m_e * constants.c**2)) +) +dEdx = n_i * constants.alpha * r_e**2 * Z**2 * (T1 + m_e_eV) * phirad print(f"dE/dx analytic = {dEdx}") -print(F"dE/dx simulated = {dEdx_simulation[-1]}") +print(f"dE/dx simulated = {dEdx_simulation[-1]}") print(f"dE/dx simulated/analytic = {dEdx_simulation[-1]/dEdx}") -assert np.abs(dEdx_simulation[-1]/dEdx - 1.) < 0.03 +assert np.abs(dEdx_simulation[-1] / dEdx - 1.0) < 0.03 sigam_total = 1.698e-28 # Calculated from table with k cutoff=1.e-4 -N_photon = n_e*n_i*L*beta*constants.c*sigam_total*dt -print(f"New photons per step simulated/analytic = {particle_number[-1,-1]/(particle_energy[-1,0]*N_photon)}") -assert np.abs(particle_number[-1,-1]/(particle_energy[-1,0]*N_photon) - 1.) < 0.03 - - +N_photon = n_e * n_i * L * beta * constants.c * sigam_total * dt +print( + f"New photons per step simulated/analytic = {particle_number[-1,-1]/(particle_energy[-1,0]*N_photon)}" +) +assert ( + np.abs(particle_number[-1, -1] / (particle_energy[-1, 0] * N_photon) - 1.0) < 0.03 +) diff --git a/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung index ba2955a5f6d..cf02199689b 100644 --- a/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung +++ b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung @@ -109,4 +109,3 @@ photon_energy.bin_min = log10(10.) photon_energy.bin_max = log10(T1) photon_energy.intervals = 10 photon_energy.precision = 18 - From ed089702317f5e729a6b0fa2def0a679e5d96303 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 13 Jan 2025 13:09:31 -0800 Subject: [PATCH 17/38] Fix missing semicolon --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 58e3ef27048..f86d354602b 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -33,7 +33,7 @@ BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, Multi bool create_photons = true; pp_collision_name.query("create_photons", create_photons); - m_exe.m_create_photons = create_photons + m_exe.m_create_photons = create_photons; amrex::ParticleReal multiplier = 1._prt; pp_collision_name.query("multiplier", multiplier); From 180a51e970895888c1e22067f9b7f298bbe9b26d Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 13 Jan 2025 13:10:48 -0800 Subject: [PATCH 18/38] Rewrite the interpolation in Photon_energy --- .../Bremsstrahlung/BremsstrahlungFunc.H | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 584763e4d80..8f39cc2e9fe 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -284,8 +284,8 @@ public: dk = (m_koT1_grid[i] - koT1_grid_im1); kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i]; amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt; - amrex::ParticleReal const next_sigma = sigma + dsigdk*dk/sigma_total; - if (random_number > next_sigma) { + amrex::ParticleReal const next_sigma = sigma + dsigdk*dk; + if (random_number*sigma_total > next_sigma) { sigma = next_sigma; koT1_grid_im1 = m_koT1_grid[i]; kdsigdk_im1 = kdsigdk_i; @@ -294,11 +294,10 @@ public: } } - amrex::ParticleReal const nnorm_i = kdsigdk_im1/koT1_grid_im1/sigma_total; - amrex::ParticleReal const nnorm_ip1 = kdsigdk_i/m_koT1_grid[i]/sigma_total; - amrex::ParticleReal const result = ((std::sqrt(nnorm_i*nnorm_i + - 2._prt*(nnorm_ip1 - nnorm_i)*(random_number - sigma)/dk) - - nnorm_i)/(nnorm_ip1 - nnorm_i))*dk + koT1_grid_im1; + amrex::ParticleReal const f_i = kdsigdk_im1/koT1_grid_im1; + amrex::ParticleReal const f_ip1 = kdsigdk_i/m_koT1_grid[i]; + amrex::ParticleReal const result = ((std::sqrt(f_i*f_i + 2._prt*(f_ip1 - f_i)*(random_number*sigma_total - sigma)/dk) + - f_i)/(f_ip1 - f_i))*dk + koT1_grid_im1; return result; } From 08088151145c84a4952aabcb42aefd299307407e Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 13 Jan 2025 13:11:26 -0800 Subject: [PATCH 19/38] Fix setting of p_mask when photons are not created --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 1 - 1 file changed, 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 8f39cc2e9fe..76c8cb08c70 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -382,7 +382,6 @@ public: // Get a random number amrex::ParticleReal const random_number = amrex::Random(engine); - p_mask[pair_index] = (random_number < q12); if (random_number < q12) { // compute weight for photon From 21eaf3183e81907c10a76c96a02169ef96153c56 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 13 Jan 2025 13:39:30 -0800 Subject: [PATCH 20/38] Another small cleanup --- .../Bremsstrahlung/BremsstrahlungFunc.H | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 76c8cb08c70..7ba9e57f514 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -294,10 +294,12 @@ public: } } - amrex::ParticleReal const f_i = kdsigdk_im1/koT1_grid_im1; - amrex::ParticleReal const f_ip1 = kdsigdk_i/m_koT1_grid[i]; - amrex::ParticleReal const result = ((std::sqrt(f_i*f_i + 2._prt*(f_ip1 - f_i)*(random_number*sigma_total - sigma)/dk) - - f_i)/(f_ip1 - f_i))*dk + koT1_grid_im1; + // k will be between k_im1 and k_i + amrex::ParticleReal const f_im1 = kdsigdk_im1/koT1_grid_im1; + amrex::ParticleReal const fi = kdsigdk_i/m_koT1_grid[i]; + amrex::ParticleReal const x = (std::sqrt(f_im1*f_im1 + 2._prt*(fi - f_im1)*(random_number*sigma_total - sigma)/dk) + - f_im1)/(fi - f_im1); + amrex::ParticleReal const result = x*dk + koT1_grid_im1; return result; } From 2548c276d3700f0d77d6ca11adbf5c27b7365430 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 10:30:52 -0800 Subject: [PATCH 21/38] Small fixes --- .../Bremsstrahlung/BremsstrahlungFunc.H | 10 +++++----- .../Bremsstrahlung/BremsstrahlungFunc.cpp | 6 +++--- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 7ba9e57f514..461f0925f19 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -179,11 +179,11 @@ public: * \param[in] E electron energy in Joules * \param[in] n1 the electron number density * \param[in] m1 the mass of the electrons - * \param[in] index the index for the electron energy grid - * \param[in] i0_cut the index below the cut off for the photon energy grid - * \param[in] koT1_cut the k of the photon energy cut off - * \param[in] kdsigdk_cut the ksigdk at the cut off - * \param[in] w0 the grid weighting at the cut off + * \param[out] index the index for the electron energy grid + * \param[out] i0_cut the index below the cut off for the photon energy grid + * \param[out] koT1_cut the k of the photon energy cut off + * \param[out] kdsigdk_cut the ksigdk at the cut off + * \param[out] w0 the grid weighting at the cut off * \result the cross section */ AMREX_GPU_HOST_DEVICE AMREX_INLINE diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index f86d354602b..49e050c5027 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -60,10 +60,10 @@ BremsstrahlungFunc::UploadCrossSection (int Z) // Convert Seltzer and Berger energy-weighted differential cross section to units of [m^2] for (int iee=0; iee < Executor::nKE; iee++) { amrex::ParticleReal const E = m_exe.m_KEgrid_eV[iee]/m_e_eV; - amrex::ParticleReal const gamma = 1.0_rt + E; + amrex::ParticleReal const gamma = 1.0_prt + E; /* betaSq = 1.0 - 1.0/gamma/gamma */ - amrex::ParticleReal const betaSq = (E*E + 2._rt*E)/gamma/gamma; - amrex::ParticleReal const scale_factor = 1.0e-31_rt*Z*Z/betaSq; + amrex::ParticleReal const betaSq = (E*E + 2._prt*E)/gamma/gamma; + amrex::ParticleReal const scale_factor = 1.0e-31_prt*Z*Z/betaSq; for (int iep=0; iep < Executor::nkoT1; iep++) { m_exe.m_kdsigdk[iee][iep] = kdsigdk[iee][iep]*scale_factor; } From 32a2a2fb46427aaa6e908945e9e97518c575103f Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 13:47:16 -0800 Subject: [PATCH 22/38] Only scale kdsigdk_cut for io_cut == 0 --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 461f0925f19..99c2b0b89df 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -237,8 +237,10 @@ public: amrex::ParticleReal const w01 = 1.0_prt - w00; kdsigdk_cut = ( w00*(w0*m_kdsigdk[index][i0_cut ] + (1.0_prt - w0)*m_kdsigdk[index+1][i0_cut]) + w01*(w0*m_kdsigdk[index][i0_cut+1] + (1.0_prt - w0)*m_kdsigdk[index+1][i0_cut+1])); - // Scale kdsigdk_cut by the centered k to mitigate divergence effect for k->0 - kdsigdk_cut *= koT1_cut/((m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt); + if (i0_cut == 0) { + // Scale kdsigdk_cut by the centered k to mitigate divergence effect for k->0 + kdsigdk_cut *= koT1_cut/((m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt); + } // Integrate the distribution using trapezoidal rule amrex::ParticleReal koT1_grid_im1 = koT1_cut; From bb265a76376e3f72a7156858b1f88c3d33c60423 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 15:04:19 -0800 Subject: [PATCH 23/38] Add comment --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 49e050c5027..6cb8ce3b93e 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -63,6 +63,7 @@ BremsstrahlungFunc::UploadCrossSection (int Z) amrex::ParticleReal const gamma = 1.0_prt + E; /* betaSq = 1.0 - 1.0/gamma/gamma */ amrex::ParticleReal const betaSq = (E*E + 2._prt*E)/gamma/gamma; + // 1.0e-31 converts mBarn to m**2 amrex::ParticleReal const scale_factor = 1.0e-31_prt*Z*Z/betaSq; for (int iep=0; iep < Executor::nkoT1; iep++) { m_exe.m_kdsigdk[iee][iep] = kdsigdk[iee][iep]*scale_factor; From f2f34fcf98e682211a6c27685860c1201065b05b Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 15:25:11 -0800 Subject: [PATCH 24/38] Fix typo --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 99c2b0b89df..d36eaceaebd 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -298,9 +298,9 @@ public: // k will be between k_im1 and k_i amrex::ParticleReal const f_im1 = kdsigdk_im1/koT1_grid_im1; - amrex::ParticleReal const fi = kdsigdk_i/m_koT1_grid[i]; - amrex::ParticleReal const x = (std::sqrt(f_im1*f_im1 + 2._prt*(fi - f_im1)*(random_number*sigma_total - sigma)/dk) - - f_im1)/(fi - f_im1); + amrex::ParticleReal const f_i = kdsigdk_i/m_koT1_grid[i]; + amrex::ParticleReal const x = (std::sqrt(f_im1*f_im1 + 2._prt*(f_i - f_im1)*(random_number*sigma_total - sigma)/dk) + - f_im1)/(f_i - f_im1); amrex::ParticleReal const result = x*dk + koT1_grid_im1; return result; From 56d504264852d3c402edc7aa9951f9c4ab48601c Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 15:31:39 -0800 Subject: [PATCH 25/38] Add precalculated cross section Uses when the calculated cut off is below the default value. --- Docs/source/usage/parameters.rst | 4 +++ .../Bremsstrahlung/BremsstrahlungFunc.H | 34 ++++++++++++------- .../Bremsstrahlung/BremsstrahlungFunc.cpp | 29 ++++++++++++++-- 3 files changed, 52 insertions(+), 15 deletions(-) diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index a91cb03c815..a5193d7f163 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -2207,6 +2207,10 @@ Details about the collision models can be found in the :ref:`theory section .create_photons`` (`integer`) Only for ``bremsstrahlung``. Whether photons will be created, defaults to 1 (true). +* ``.koT1_cut`` (`float`) + Only for ``bremsstrahlung``. Minimum energy of the photons created. + This is relative to the electron energy, defaulting to 1.e-4. + .. _running-cpp-parameters-numerics: Numerics and algorithms diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index d36eaceaebd..bd9ac7fd38b 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -207,7 +207,7 @@ public: } // find lo-index for koT1 cutoff (will typically be 1) - koT1_cut = std::max(KEcut_eV/KErel_eV, 1.0e-4_prt); + koT1_cut = std::max(KEcut_eV/KErel_eV, m_koT1_cut_default); i0_cut = 0; for (int i=1; i < nkoT1; i++) { if (m_koT1_grid[i] > koT1_cut) { @@ -242,19 +242,25 @@ public: kdsigdk_cut *= koT1_cut/((m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt); } - // Integrate the distribution using trapezoidal rule - amrex::ParticleReal koT1_grid_im1 = koT1_cut; - amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; - amrex::ParticleReal sigma = 0.0_prt; - for (int i=i0_cut+1; i < nkoT1; i++) { - amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1); - amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i]; - amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt; - sigma = sigma + dsigdk*dk; - koT1_grid_im1 = m_koT1_grid[i]; - kdsigdk_im1 = kdsigdk_i; + amrex::ParticleReal sigma_total; + if (KEcut_eV/KErel_eV <= m_koT1_cut_default) { + // Use precalculated value + sigma_total = w0*m_sigma_total[index] + (1.0_prt - w0)*m_sigma_total[index+1]; + } else { + // Integrate the distribution using trapezoidal rule + amrex::ParticleReal koT1_grid_im1 = koT1_cut; + amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; + amrex::ParticleReal sigma = 0.0_prt; + for (int i=i0_cut+1; i < nkoT1; i++) { + amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1); + amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i]; + amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt; + sigma = sigma + dsigdk*dk; + koT1_grid_im1 = m_koT1_grid[i]; + kdsigdk_im1 = kdsigdk_i; + } + sigma_total = sigma; } - amrex::ParticleReal const sigma_total = sigma; // total cross section [m^2] return sigma_total; } @@ -443,6 +449,7 @@ public: bool m_create_photons = true; bool m_need_product_data = true; amrex::ParticleReal m_multiplier; + amrex::ParticleReal m_koT1_cut_default; // relative photon energy grid ( koT1 = Ephoton_eV/KErel_eV ) amrex::GpuArray m_koT1_grid = {0., 0.025, 0.05, 0.075, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.97, 0.99, 1.00}; @@ -452,6 +459,7 @@ public: // differential total cross section amrex::GpuArray, nKE> m_kdsigdk; + amrex::GpuArray m_sigma_total; }; diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 6cb8ce3b93e..4eff0bcf84f 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -41,6 +41,10 @@ BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, Multi "BremsstrahlungFunc: The multiplier must be greater than or equal to one"); m_exe.m_multiplier = multiplier; + amrex::ParticleReal koT1_cut_default = 1.e-4; + pp_collision_name.query("koT1_cut", koT1_cut_default); + m_exe.m_koT1_cut_default = koT1_cut_default; + // Fill in the m_kdsigdk array UploadCrossSection(Z); @@ -57,6 +61,11 @@ BremsstrahlungFunc::UploadCrossSection (int Z) std::vector> & kdsigdk = m_kdsigdk_map.at(Z); + // kdsigdk at the default cutoff + int i0_cut = 0; + amrex::ParticleReal const koT1_cut = m_exe.m_koT1_cut_default; + amrex::ParticleReal const w00 = (m_exe.m_koT1_grid[i0_cut+1] - koT1_cut)/(m_exe.m_koT1_grid[i0_cut+1] - m_exe.m_koT1_grid[i0_cut]); + // Convert Seltzer and Berger energy-weighted differential cross section to units of [m^2] for (int iee=0; iee < Executor::nKE; iee++) { amrex::ParticleReal const E = m_exe.m_KEgrid_eV[iee]/m_e_eV; @@ -65,9 +74,25 @@ BremsstrahlungFunc::UploadCrossSection (int Z) amrex::ParticleReal const betaSq = (E*E + 2._prt*E)/gamma/gamma; // 1.0e-31 converts mBarn to m**2 amrex::ParticleReal const scale_factor = 1.0e-31_prt*Z*Z/betaSq; - for (int iep=0; iep < Executor::nkoT1; iep++) { - m_exe.m_kdsigdk[iee][iep] = kdsigdk[iee][iep]*scale_factor; + for (int i=0; i < Executor::nkoT1; i++) { + m_exe.m_kdsigdk[iee][i] = kdsigdk[iee][i]*scale_factor; + } + + // Calculate the total cross section using the default k-cutoff + amrex::ParticleReal const kdsigdk_cut = w00*m_exe.m_kdsigdk[iee][i0_cut] + (1.0_prt - w00)*m_exe.m_kdsigdk[iee][i0_cut+1]; + amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut*koT1_cut/((m_exe.m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt); + amrex::ParticleReal koT1_im1 = koT1_cut; + amrex::ParticleReal sigma = 0.0_prt; + for (int i=i0_cut+1; i < Executor::nkoT1; i++) { + amrex::ParticleReal const koT1_i = m_exe.m_koT1_grid[i]; + amrex::ParticleReal const dk = (koT1_i - koT1_im1); + amrex::ParticleReal const kdsigdk_i = m_exe.m_kdsigdk[iee][i]; + amrex::ParticleReal const dsigdk = (kdsigdk_i/koT1_i + kdsigdk_im1/koT1_im1)*0.5_prt; + sigma = sigma + dsigdk*dk; + koT1_im1 = koT1_i; + kdsigdk_im1 = kdsigdk_i; } + m_exe.m_sigma_total[iee] = sigma; } } From 98bc9e8eab784e076671ae28a55d67be253a5cfe Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 14 Jan 2025 23:33:36 +0000 Subject: [PATCH 26/38] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../Tests/collision/analysis_collision_1d_Bremsstrahlung.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py index d3242d6c91e..508f8fbff6a 100755 --- a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py +++ b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py @@ -94,13 +94,13 @@ print(f"dE/dx analytic = {dEdx}") print(f"dE/dx simulated = {dEdx_simulation[-1]}") -print(f"dE/dx simulated/analytic = {dEdx_simulation[-1]/dEdx}") +print(f"dE/dx simulated/analytic = {dEdx_simulation[-1] / dEdx}") assert np.abs(dEdx_simulation[-1] / dEdx - 1.0) < 0.03 sigam_total = 1.698e-28 # Calculated from table with k cutoff=1.e-4 N_photon = n_e * n_i * L * beta * constants.c * sigam_total * dt print( - f"New photons per step simulated/analytic = {particle_number[-1,-1]/(particle_energy[-1,0]*N_photon)}" + f"New photons per step simulated/analytic = {particle_number[-1, -1] / (particle_energy[-1, 0] * N_photon)}" ) assert ( np.abs(particle_number[-1, -1] / (particle_energy[-1, 0] * N_photon) - 1.0) < 0.03 From eeadf9a3e8d99552846bf53e1ef5f955f2d95504 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 16:19:13 -0800 Subject: [PATCH 27/38] Clean up in analysis_collision_1d_Bremsstrahlung.py --- .../Tests/collision/analysis_collision_1d_Bremsstrahlung.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py index 508f8fbff6a..d8a3ca565ab 100755 --- a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py +++ b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py @@ -29,9 +29,6 @@ particle_number = np.loadtxt( os.path.join("diags", "reducedfiles", "particle_number.txt"), skiprows=1 ) -photon_energy = np.loadtxt( - os.path.join("diags", "reducedfiles", "photon_energy.txt"), skiprows=1 -) total_energy = particle_energy[:, 2] electron_energy = particle_energy[:, 3] From 54ff0532c58ca7a790d90ab77c073e64d4048df7 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 14 Jan 2025 16:19:55 -0800 Subject: [PATCH 28/38] Clean up inputs_test_1d_collision_z_Bremsstrahlung --- .../inputs_test_1d_collision_z_Bremsstrahlung | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung index cf02199689b..dd6156c6510 100644 --- a/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung +++ b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung @@ -94,18 +94,10 @@ diagnostics.diags_names = diag1 diag1.intervals = 100 diag1.diag_type = Full -warpx.reduced_diags_names = particle_energy particle_number photon_energy +warpx.reduced_diags_names = particle_energy particle_number particle_energy.type = ParticleEnergy particle_energy.intervals = 10 particle_energy.precision = 18 particle_number.type = ParticleNumber particle_number.intervals = 10 particle_number.precision = 18 -photon_energy.type = ParticleHistogram -photon_energy.species = photons -photon_energy.histogram_function(t,x,y,z,ux,uy,uz) = log10(sqrt(ux*ux+uy*uy+uz*uz)*clight*clight/q_e*m_e) -photon_energy.bin_number = 100 -photon_energy.bin_min = log10(10.) -photon_energy.bin_max = log10(T1) -photon_energy.intervals = 10 -photon_energy.precision = 18 From f9756fdc850b617cb5921d658c4eabe3e3ae5557 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 20 Jan 2025 08:58:30 -0800 Subject: [PATCH 29/38] Revert to using the midpoint value of k This gives results that more closely matches the phirad from Seltzer and Berger --- .../test_1d_collision_z_Bremsstrahlung.json | 32 +++++++++---------- .../Bremsstrahlung/BremsstrahlungFunc.H | 17 +++++----- .../Bremsstrahlung/BremsstrahlungFunc.cpp | 7 ++-- 3 files changed, 28 insertions(+), 28 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json b/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json index d8edba3a5b4..6adbbcc826f 100644 --- a/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json +++ b/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json @@ -10,25 +10,25 @@ "jy": 0.0, "jz": 0.0 }, + "ions": { + "particle_momentum_x": 6.244058619001046e-18, + "particle_momentum_y": 6.266839027572808e-18, + "particle_momentum_z": 6.237935835816471e-18, + "particle_position_x": 0.0511999988707094, + "particle_weight": 5.470000000000001e+25 + }, "electrons": { - "particle_momentum_x": 1.7785184549415273e-26, - "particle_momentum_y": 1.769780531165544e-26, - "particle_momentum_z": 7.780435420929519e-17, - "particle_position_x": 0.0511996620835328, + "particle_momentum_x": 1.8547285435878255e-26, + "particle_momentum_y": 1.8336044974219876e-26, + "particle_momentum_z": 7.780403913830049e-17, + "particle_position_x": 0.05119965192677569, "particle_weight": 5.470000000000001e+24 }, "photons": { - "particle_momentum_x": 3.9236321845123015e-27, - "particle_momentum_y": 3.892475715179289e-27, - "particle_momentum_z": 1.2828655290730977e-18, - "particle_position_x": 0.013297153040958081, - "particle_weight": 1.4269541992187505e+22 - }, - "ions": { - "particle_momentum_x": 6.245889688274839e-18, - "particle_momentum_y": 6.2336526718397175e-18, - "particle_momentum_z": 6.258023561029987e-18, - "particle_position_x": 0.05120000226548876, - "particle_weight": 5.470000000000001e+25 + "particle_momentum_x": 4.200876982473409e-27, + "particle_momentum_y": 4.172037413327287e-27, + "particle_momentum_z": 1.3125096507793955e-18, + "particle_position_x": 0.013402699474072983, + "particle_weight": 1.4387061523437505e+22 } } diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index bd9ac7fd38b..6af95c1f986 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -237,10 +237,6 @@ public: amrex::ParticleReal const w01 = 1.0_prt - w00; kdsigdk_cut = ( w00*(w0*m_kdsigdk[index][i0_cut ] + (1.0_prt - w0)*m_kdsigdk[index+1][i0_cut]) + w01*(w0*m_kdsigdk[index][i0_cut+1] + (1.0_prt - w0)*m_kdsigdk[index+1][i0_cut+1])); - if (i0_cut == 0) { - // Scale kdsigdk_cut by the centered k to mitigate divergence effect for k->0 - kdsigdk_cut *= koT1_cut/((m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt); - } amrex::ParticleReal sigma_total; if (KEcut_eV/KErel_eV <= m_koT1_cut_default) { @@ -252,9 +248,10 @@ public: amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; amrex::ParticleReal sigma = 0.0_prt; for (int i=i0_cut+1; i < nkoT1; i++) { - amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1); amrex::ParticleReal const kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i]; - amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt; + amrex::ParticleReal const dk = (m_koT1_grid[i] - koT1_grid_im1); + amrex::ParticleReal const k_ave = (m_koT1_grid[i] + koT1_grid_im1)*0.5_prt; + amrex::ParticleReal const dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/k_ave; sigma = sigma + dsigdk*dk; koT1_grid_im1 = m_koT1_grid[i]; kdsigdk_im1 = kdsigdk_i; @@ -287,11 +284,13 @@ public: amrex::ParticleReal sigma = 0._prt; amrex::ParticleReal kdsigdk_i; amrex::ParticleReal dk; + amrex::ParticleReal k_ave; int i; for (i=i0_cut+1; i < nkoT1; i++) { dk = (m_koT1_grid[i] - koT1_grid_im1); + k_ave = (m_koT1_grid[i] + koT1_grid_im1)*0.5_prt; kdsigdk_i = w0*m_kdsigdk[index][i] + (1.0_prt - w0)*m_kdsigdk[index+1][i]; - amrex::ParticleReal const dsigdk = (kdsigdk_i/m_koT1_grid[i] + kdsigdk_im1/koT1_grid_im1)*0.5_prt; + amrex::ParticleReal const dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/k_ave; amrex::ParticleReal const next_sigma = sigma + dsigdk*dk; if (random_number*sigma_total > next_sigma) { sigma = next_sigma; @@ -303,8 +302,8 @@ public: } // k will be between k_im1 and k_i - amrex::ParticleReal const f_im1 = kdsigdk_im1/koT1_grid_im1; - amrex::ParticleReal const f_i = kdsigdk_i/m_koT1_grid[i]; + amrex::ParticleReal const f_im1 = kdsigdk_im1/k_ave; + amrex::ParticleReal const f_i = kdsigdk_i/k_ave; amrex::ParticleReal const x = (std::sqrt(f_im1*f_im1 + 2._prt*(f_i - f_im1)*(random_number*sigma_total - sigma)/dk) - f_im1)/(f_i - f_im1); amrex::ParticleReal const result = x*dk + koT1_grid_im1; diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 4eff0bcf84f..0dff7bffbeb 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -80,14 +80,15 @@ BremsstrahlungFunc::UploadCrossSection (int Z) // Calculate the total cross section using the default k-cutoff amrex::ParticleReal const kdsigdk_cut = w00*m_exe.m_kdsigdk[iee][i0_cut] + (1.0_prt - w00)*m_exe.m_kdsigdk[iee][i0_cut+1]; - amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut*koT1_cut/((m_exe.m_koT1_grid[i0_cut+1] + koT1_cut)*0.5_prt); + amrex::ParticleReal kdsigdk_im1 = kdsigdk_cut; amrex::ParticleReal koT1_im1 = koT1_cut; amrex::ParticleReal sigma = 0.0_prt; for (int i=i0_cut+1; i < Executor::nkoT1; i++) { amrex::ParticleReal const koT1_i = m_exe.m_koT1_grid[i]; - amrex::ParticleReal const dk = (koT1_i - koT1_im1); amrex::ParticleReal const kdsigdk_i = m_exe.m_kdsigdk[iee][i]; - amrex::ParticleReal const dsigdk = (kdsigdk_i/koT1_i + kdsigdk_im1/koT1_im1)*0.5_prt; + amrex::ParticleReal const dk = (koT1_i - koT1_im1); + amrex::ParticleReal const k_ave = (koT1_i + koT1_im1)*0.5_prt; + amrex::ParticleReal const dsigdk = (kdsigdk_i + kdsigdk_im1)*0.5_prt/k_ave; sigma = sigma + dsigdk*dk; koT1_im1 = koT1_i; kdsigdk_im1 = kdsigdk_i; From 7a91252bbffc123b735a5b3a1d882ab67e0fcc8f Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 20 Jan 2025 09:00:01 -0800 Subject: [PATCH 30/38] Fix the sign of the photon energy --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 6af95c1f986..db92f1c6056 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -427,8 +427,8 @@ public: amrex::ParticleReal const KE_after = m_e_eV*gbp1sq_prime/(1.0_prt + gammap1_prime); // update energy lost to bremsstrahlung - amrex::ParticleReal const deltaE_eV = KE_after - KE_before; - /* m_deltaE_bremsstrahlung -= deltaE_eV*weight1*PhysConst::q_e; // Joules */ + amrex::ParticleReal const deltaE_eV = KE_before - KE_after; + /* m_deltaE_bremsstrahlung += deltaE_eV*weight1*PhysConst::q_e; // Joules */ if (m_create_photons) { // compute photon energy in lab frame From 2b087e831beea481868784ac815bf411b3a6f255 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 20 Jan 2025 11:06:06 -0800 Subject: [PATCH 31/38] Fix const --- .../BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 0dff7bffbeb..6ac5d4a4017 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -62,7 +62,7 @@ BremsstrahlungFunc::UploadCrossSection (int Z) std::vector> & kdsigdk = m_kdsigdk_map.at(Z); // kdsigdk at the default cutoff - int i0_cut = 0; + int const i0_cut = 0; amrex::ParticleReal const koT1_cut = m_exe.m_koT1_cut_default; amrex::ParticleReal const w00 = (m_exe.m_koT1_grid[i0_cut+1] - koT1_cut)/(m_exe.m_koT1_grid[i0_cut+1] - m_exe.m_koT1_grid[i0_cut]); From 326734494832de7161f20fc80df17ef780d80893 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 21 Jan 2025 09:38:15 -0800 Subject: [PATCH 32/38] Fix CI analysis script --- .../Tests/collision/analysis_collision_1d_Bremsstrahlung.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py index d8a3ca565ab..bca6ed3d75f 100755 --- a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py +++ b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py @@ -94,11 +94,11 @@ print(f"dE/dx simulated/analytic = {dEdx_simulation[-1] / dEdx}") assert np.abs(dEdx_simulation[-1] / dEdx - 1.0) < 0.03 -sigam_total = 1.698e-28 # Calculated from table with k cutoff=1.e-4 -N_photon = n_e * n_i * L * beta * constants.c * sigam_total * dt +sigma_total = 1.818e-28 # Calculated from table with k cutoff=1.e-4 +N_photon = n_e * n_i * L * beta * constants.c * sigma_total * dt print( f"New photons per step simulated/analytic = {particle_number[-1, -1] / (particle_energy[-1, 0] * N_photon)}" ) assert ( - np.abs(particle_number[-1, -1] / (particle_energy[-1, 0] * N_photon) - 1.0) < 0.03 + np.abs(particle_number[-1, -1] / (particle_energy[-1, 0] * N_photon) - 1.0) < 0.02 ) From d9b05ca8b2e03df36ddeb5cb7749eae21649bc84 Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 21 Jan 2025 15:05:58 -0800 Subject: [PATCH 33/38] Clean up auto in Bremsstrahlung/PhotonCreationFunc --- .../BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H index 77e2bb42b1e..920460095e8 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H @@ -163,13 +163,13 @@ public: // Set the weight of the new particle to p_pair_reaction_weight[i] soa_products_data[0].m_rdata[PIdx::w][product_index] = p_pair_reaction_weight[i]; - auto& ux1 = soa_products_data[0].m_rdata[PIdx::ux][product_index]; - auto& uy1 = soa_products_data[0].m_rdata[PIdx::uy][product_index]; - auto& uz1 = soa_products_data[0].m_rdata[PIdx::uz][product_index]; + amrex::ParticleReal & ux1 = soa_products_data[0].m_rdata[PIdx::ux][product_index]; + amrex::ParticleReal & uy1 = soa_products_data[0].m_rdata[PIdx::uy][product_index]; + amrex::ParticleReal & uz1 = soa_products_data[0].m_rdata[PIdx::uz][product_index]; // Normalize out the electron velocity and multiply by the photon momentum, E/c // Also, the photon momentum is normalized by m_e - auto u1 = std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1); + amrex::ParticleReal const u1 = std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1); ux1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; uy1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; uz1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; From 0f3e8817a0466095b5b41a95fb894487c09f125a Mon Sep 17 00:00:00 2001 From: David Grote Date: Tue, 21 Jan 2025 16:02:49 -0800 Subject: [PATCH 34/38] Ions are updated, now conserving energy and momentum --- .../Bremsstrahlung/BremsstrahlungFunc.H | 100 ++++++++++-------- 1 file changed, 58 insertions(+), 42 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index db92f1c6056..29e77d4ddb5 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -90,7 +90,7 @@ public: amrex::ParticleReal const n1, amrex::ParticleReal const n2, amrex::ParticleReal const /*T1*/, amrex::ParticleReal const /*T2*/, amrex::ParticleReal const /*q1*/, amrex::ParticleReal const /*q2*/, - amrex::ParticleReal const m1, amrex::ParticleReal const /*m2*/, + amrex::ParticleReal const m1, amrex::ParticleReal const m2, amrex::Real const dt, amrex::Real const /*dV*/, index_type coll_idx, index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask, index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2, @@ -150,7 +150,7 @@ public: BremsstrahlungEvent( u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], - m1, weight1[ I1[i1] ], weight2[ I2[i2] ], + m1, m2, weight1[ I1[i1] ], weight2[ I2[i2] ], n1, n2, dt, static_cast(pair_index), p_mask, p_pair_reaction_weight, p_product_data, engine); @@ -315,7 +315,7 @@ public: * * \param[in] u1x, u1y, u1z the gamma*velocity of the first species, the electrons * \param[in] u2x, u2y, u2z the gamma*velocity of the second species, the atoms - * \param[in] m1 the mass of the electrons + * \param[in] m1, m2 the mass of the electrons and ions * \param[in] weight1 the particle weight of the electrons * \param[in] weight2 the particle weight of the target (unused) * \param[in] n1 the number density of the electrons @@ -329,10 +329,10 @@ public: */ AMREX_GPU_HOST_DEVICE AMREX_INLINE void BremsstrahlungEvent (amrex::ParticleReal & u1x, amrex::ParticleReal & u1y, - amrex::ParticleReal & u1z, amrex::ParticleReal const& u2x, - amrex::ParticleReal const& u2y, amrex::ParticleReal const& u2z, - amrex::ParticleReal m1, - amrex::ParticleReal weight1, amrex::ParticleReal /*weight2*/, + amrex::ParticleReal & u1z, amrex::ParticleReal & u2x, + amrex::ParticleReal & u2y, amrex::ParticleReal & u2z, + amrex::ParticleReal m1, amrex::ParticleReal m2, + amrex::ParticleReal weight1, amrex::ParticleReal weight2, amrex::ParticleReal n1, amrex::ParticleReal n2, amrex::Real dt, int pair_index, index_type* AMREX_RESTRICT p_mask, @@ -343,13 +343,12 @@ public: using namespace amrex::literals; constexpr auto c2 = PhysConst::c * PhysConst::c; - amrex::ParticleReal const m_e_eV = m1*PhysConst::c*PhysConst::c/PhysConst::q_e; // 0.511e6 // compute gamma for particles 1 (electron) and 2 (ion) - amrex::ParticleReal const gbp1sq = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; - amrex::ParticleReal const gbp2sq = (u2x*u2x + u2y*u2y + u2z*u2z)/c2; - amrex::ParticleReal const gammap1 = std::sqrt(1.0_prt + gbp1sq); - amrex::ParticleReal const gammap2 = std::sqrt(1.0_prt + gbp2sq); + amrex::ParticleReal const u1sq = (u1x*u1x + u1y*u1y + u1z*u1z); + amrex::ParticleReal const u2sq = (u2x*u2x + u2y*u2y + u2z*u2z); + amrex::ParticleReal const gamma1 = std::sqrt(1.0_prt + u1sq/c2); + amrex::ParticleReal const gamma2 = std::sqrt(1.0_prt + u2sq/c2); // transform electron to frame where ion is at rest amrex::ParticleReal u1x_rel = u1x; @@ -358,9 +357,9 @@ public: ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, u2x, u2y, u2z); // compute electron KE in frame where ion is at rest - amrex::ParticleReal const gbp1sq_rel = (u1x_rel*u1x_rel + u1y_rel*u1y_rel + u1z_rel*u1z_rel)/c2; - amrex::ParticleReal const gammap1_rel = std::sqrt(1.0_prt + gbp1sq_rel); - amrex::ParticleReal const KE_eV = m_e_eV*gbp1sq_rel/(1.0_prt + gammap1_rel); + amrex::ParticleReal const u1sq_rel = (u1x_rel*u1x_rel + u1y_rel*u1y_rel + u1z_rel*u1z_rel); + amrex::ParticleReal const gamma1_rel = std::sqrt(1.0_prt + u1sq_rel/c2); + amrex::ParticleReal const KE_eV = m1*u1sq_rel/(1.0_prt + gamma1_rel)/PhysConst::q_e; int i0_cut, index; amrex::ParticleReal koT1_cut, kdsigdk_cut; @@ -372,13 +371,13 @@ public: // determine if the pair collide and then do the collision amrex::ParticleReal fmulti = m_multiplier; // >= 1.0 - amrex::ParticleReal const gamma_factor = gammap1_rel/(gammap1*gammap2); - amrex::ParticleReal const beta1 = std::sqrt(gbp1sq_rel)/gammap1_rel; // magnitude electron beta - amrex::ParticleReal arg = fmulti*beta1*PhysConst::c*sigma_total*n2*dt*gamma_factor; + amrex::ParticleReal const gamma_factor = gamma1_rel/(gamma1*gamma2); + amrex::ParticleReal const v1 = std::sqrt(u1sq_rel)/gamma1_rel; // magnitude electron velocity + amrex::ParticleReal arg = fmulti*v1*sigma_total*n2*dt*gamma_factor; if (arg > 1.0_prt) { #ifndef AMREX_USE_GPU amrex::Print() << "Notice: arg = " << arg << " in interSpeciesBremsstrahlung" << "\n"; - amrex::Print() << " beta1 = " << beta1 << "\n"; + amrex::Print() << " v1 = " << v1 << "\n"; amrex::Print() << " n2 = " << n2 << "\n"; amrex::Print() << " sigma_total = " << sigma_total << "\n"; #endif @@ -396,43 +395,60 @@ public: // compute weight for photon amrex::ParticleReal const w_photon = weight1/fmulti; - // compute lab-frame electron kinetic energy before emission - amrex::ParticleReal const KE_before = m_e_eV*gbp1sq/(1.0_prt + gammap1); - // get energy of photon (hbar*omega) amrex::ParticleReal const random_number2 = amrex::Random(engine); amrex::ParticleReal const EphotonoT1 = Photon_energy(index, i0_cut, koT1_cut, kdsigdk_cut, w0, sigma_total, random_number2); - amrex::ParticleReal const Ephoton_eV = EphotonoT1*KE_eV; - - // adjust electron gamma and beta - amrex::ParticleReal const KE_prime_eV = KE_eV - Ephoton_eV*w_photon/weight1; - amrex::ParticleReal const gamma_prime = 1.0_prt + KE_prime_eV/m_e_eV; - - // rescale electron beta in frame where ion is at rest - amrex::ParticleReal const u_rel_ratio = std::sqrt((gamma_prime + 1.0_prt)*KE_prime_eV/m_e_eV/gbp1sq_rel); - u1x_rel *= u_rel_ratio; - u1y_rel *= u_rel_ratio; - u1z_rel *= u_rel_ratio; - - // Lorentz transform electron back to lab frame + amrex::ParticleReal const Ephoton = EphotonoT1*KE_eV*PhysConst::q_e; + + // Calculate electron and ion energy and momentum using conservation + amrex::ParticleReal const u1_rel = std::sqrt(u1sq_rel); + amrex::ParticleReal const A = u1_rel/PhysConst::c - Ephoton*w_photon/(weight1*m1*c2); + amrex::ParticleReal const B = gamma1_rel + m2/m1 - Ephoton*w_photon/(weight1*m1*c2); + amrex::ParticleReal const D = m2*m2/(m1*m1) + A*A - B*B - 1.0_prt; + amrex::ParticleReal const a = 4.0_prt*(B*B - A*A); + amrex::ParticleReal const b = 4.0_prt*A*D; + amrex::ParticleReal const c = 4.0_prt*B*B - D*D; + // Use the "plus" solution since it keeps the electon motion forward + amrex::ParticleReal const u1_rel_after = (-b + std::sqrt(b*b - 4.0_prt*a*c))/(2.0_prt*a)*PhysConst::c; + amrex::ParticleReal const u2_rel_after = (A - u1_rel_after/PhysConst::c)*m1/m2*PhysConst::c; + + u1x_rel *= u1_rel_after/u1_rel; + u1y_rel *= u1_rel_after/u1_rel; + u1z_rel *= u1_rel_after/u1_rel; + amrex::ParticleReal u2x_rel = u2_rel_after*u1x_rel/u1_rel_after; + amrex::ParticleReal u2y_rel = u2_rel_after*u1y_rel/u1_rel_after; + amrex::ParticleReal u2z_rel = u2_rel_after*u1z_rel/u1_rel_after; + + // Lorentz transform electron and ion back to lab frame ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, -u2x, -u2y, -u2z); + ParticleUtils::doLorentzTransformWithU(u2x_rel, u2y_rel, u2z_rel, -u2x, -u2y, -u2z); u1x = u1x_rel; u1y = u1y_rel; u1z = u1z_rel; + u2x = u2x_rel; + u2y = u2y_rel; + u2z = u2z_rel; + + // compute lab-frame electron kinetic energy before emission + amrex::ParticleReal const KE1_before = weight1*m1*u1sq/(1.0_prt + gamma1); + amrex::ParticleReal const KE2_before = weight2*m2*u2sq/(1.0_prt + gamma2); - // compute electron kinetic energy after emission - amrex::ParticleReal const gbp1sq_prime = (u1x*u1x + u1y*u1y + u1z*u1z)/c2; - amrex::ParticleReal const gammap1_prime = std::sqrt(1.0_prt + gbp1sq_prime); - amrex::ParticleReal const KE_after = m_e_eV*gbp1sq_prime/(1.0_prt + gammap1_prime); + // compute electron and ion kinetic energy after emission + amrex::ParticleReal const u1sq_after = (u1x*u1x + u1y*u1y + u1z*u1z); + amrex::ParticleReal const u2sq_after = (u2x*u2x + u2y*u2y + u2z*u2z); + amrex::ParticleReal const gamma1_after = std::sqrt(1.0_prt + u1sq_after/c2); + amrex::ParticleReal const gamma2_after = std::sqrt(1.0_prt + u2sq_after/c2); + amrex::ParticleReal const KE1_after = weight1*m1*u1sq_after/(1.0_prt + gamma1_after); + amrex::ParticleReal const KE2_after = weight2*m2*u2sq_after/(1.0_prt + gamma2_after); // update energy lost to bremsstrahlung - amrex::ParticleReal const deltaE_eV = KE_before - KE_after; - /* m_deltaE_bremsstrahlung += deltaE_eV*weight1*PhysConst::q_e; // Joules */ + amrex::ParticleReal const deltaE = KE1_before + KE2_before - KE1_after - KE2_after; + /* m_deltaE_bremsstrahlung += deltaE*weight1; // Joules */ if (m_create_photons) { // compute photon energy in lab frame - amrex::ParticleReal const Ephoton_lab = deltaE_eV*weight1/w_photon*PhysConst::q_e; // Joules + amrex::ParticleReal const Ephoton_lab = deltaE/w_photon; // Joules p_mask[pair_index] = true; p_pair_reaction_weight[pair_index] = w_photon; From 0e524a3fd7b5f0bad375fb286b9808e49a625855 Mon Sep 17 00:00:00 2001 From: David Grote Date: Fri, 24 Jan 2025 16:32:33 -0800 Subject: [PATCH 35/38] Rework the collision to conserve energy and momentum Now the ion momentum is also updated --- .../analysis_collision_1d_Bremsstrahlung.py | 25 ++-- .../inputs_test_1d_collision_z_Bremsstrahlung | 5 +- .../Bremsstrahlung/BremsstrahlungFunc.H | 64 ++-------- .../Bremsstrahlung/BremsstrahlungFunc.cpp | 4 - .../Bremsstrahlung/PhotonCreationFunc.H | 111 +++++++++++++++--- .../Bremsstrahlung/PhotonCreationFunc.cpp | 6 + Source/Utils/ParticleUtils.H | 46 +++++++- 7 files changed, 167 insertions(+), 94 deletions(-) diff --git a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py index bca6ed3d75f..03f13627c21 100755 --- a/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py +++ b/Examples/Tests/collision/analysis_collision_1d_Bremsstrahlung.py @@ -26,6 +26,9 @@ particle_energy = np.loadtxt( os.path.join("diags", "reducedfiles", "particle_energy.txt"), skiprows=1 ) +particle_momentum = np.loadtxt( + os.path.join("diags", "reducedfiles", "particle_momentum.txt"), skiprows=1 +) particle_number = np.loadtxt( os.path.join("diags", "reducedfiles", "particle_number.txt"), skiprows=1 ) @@ -35,7 +38,12 @@ ion_energy = particle_energy[:, 4] photon_energy = particle_energy[:, 5] -energy_tolerance = 1.0e-12 +total_momentum = particle_momentum[:, 2] +electron_momentum = particle_momentum[:, 3] +ion_momentum = particle_momentum[:, 4] +photon_momentum = particle_momentum[:, 5] + +energy_tolerance = 1.0e-11 print(f"initial total energy = {total_energy[0]}") print(f"final total energy = {total_energy[-1]}") @@ -44,17 +52,14 @@ print(f"change in total energy = {dE_total}") assert dE_total < energy_tolerance -print(f"initial electron energy = {electron_energy[0]}") -print(f"final electron energy = {electron_energy[-1]}") -print(f"final photon energy = {photon_energy[-1]}") +momentum_tolerance = 1.0e-12 -dE_electron_photon = ( - np.abs(electron_energy[-1] + photon_energy[-1] - electron_energy[0]) - / electron_energy[0] -) -print(f"electron, photon energy change = {dE_electron_photon}") -assert dE_electron_photon < energy_tolerance +print(f"initial total momentum = {total_momentum[0]}") +print(f"final total momentum = {total_momentum[-1]}") +dP_total = np.abs(total_momentum[-1] - total_momentum[0]) / total_momentum[0] +print(f"change in total momentum = {dP_total}") +assert dP_total < momentum_tolerance print() diff --git a/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung index dd6156c6510..84f06f288cd 100644 --- a/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung +++ b/Examples/Tests/collision/inputs_test_1d_collision_z_Bremsstrahlung @@ -94,10 +94,13 @@ diagnostics.diags_names = diag1 diag1.intervals = 100 diag1.diag_type = Full -warpx.reduced_diags_names = particle_energy particle_number +warpx.reduced_diags_names = particle_energy particle_momentum particle_number particle_energy.type = ParticleEnergy particle_energy.intervals = 10 particle_energy.precision = 18 +particle_momentum.type = ParticleMomentum +particle_momentum.intervals = 10 +particle_momentum.precision = 18 particle_number.type = ParticleNumber particle_number.intervals = 10 particle_number.precision = 18 diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H index 29e77d4ddb5..9b293523ec9 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H @@ -331,8 +331,8 @@ public: void BremsstrahlungEvent (amrex::ParticleReal & u1x, amrex::ParticleReal & u1y, amrex::ParticleReal & u1z, amrex::ParticleReal & u2x, amrex::ParticleReal & u2y, amrex::ParticleReal & u2z, - amrex::ParticleReal m1, amrex::ParticleReal m2, - amrex::ParticleReal weight1, amrex::ParticleReal weight2, + amrex::ParticleReal m1, amrex::ParticleReal /*m2*/, + amrex::ParticleReal weight1, amrex::ParticleReal /*weight2*/, amrex::ParticleReal n1, amrex::ParticleReal n2, amrex::Real dt, int pair_index, index_type* AMREX_RESTRICT p_mask, @@ -401,59 +401,12 @@ public: w0, sigma_total, random_number2); amrex::ParticleReal const Ephoton = EphotonoT1*KE_eV*PhysConst::q_e; - // Calculate electron and ion energy and momentum using conservation - amrex::ParticleReal const u1_rel = std::sqrt(u1sq_rel); - amrex::ParticleReal const A = u1_rel/PhysConst::c - Ephoton*w_photon/(weight1*m1*c2); - amrex::ParticleReal const B = gamma1_rel + m2/m1 - Ephoton*w_photon/(weight1*m1*c2); - amrex::ParticleReal const D = m2*m2/(m1*m1) + A*A - B*B - 1.0_prt; - amrex::ParticleReal const a = 4.0_prt*(B*B - A*A); - amrex::ParticleReal const b = 4.0_prt*A*D; - amrex::ParticleReal const c = 4.0_prt*B*B - D*D; - // Use the "plus" solution since it keeps the electon motion forward - amrex::ParticleReal const u1_rel_after = (-b + std::sqrt(b*b - 4.0_prt*a*c))/(2.0_prt*a)*PhysConst::c; - amrex::ParticleReal const u2_rel_after = (A - u1_rel_after/PhysConst::c)*m1/m2*PhysConst::c; - - u1x_rel *= u1_rel_after/u1_rel; - u1y_rel *= u1_rel_after/u1_rel; - u1z_rel *= u1_rel_after/u1_rel; - amrex::ParticleReal u2x_rel = u2_rel_after*u1x_rel/u1_rel_after; - amrex::ParticleReal u2y_rel = u2_rel_after*u1y_rel/u1_rel_after; - amrex::ParticleReal u2z_rel = u2_rel_after*u1z_rel/u1_rel_after; - - // Lorentz transform electron and ion back to lab frame - ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, -u2x, -u2y, -u2z); - ParticleUtils::doLorentzTransformWithU(u2x_rel, u2y_rel, u2z_rel, -u2x, -u2y, -u2z); - u1x = u1x_rel; - u1y = u1y_rel; - u1z = u1z_rel; - u2x = u2x_rel; - u2y = u2y_rel; - u2z = u2z_rel; - - // compute lab-frame electron kinetic energy before emission - amrex::ParticleReal const KE1_before = weight1*m1*u1sq/(1.0_prt + gamma1); - amrex::ParticleReal const KE2_before = weight2*m2*u2sq/(1.0_prt + gamma2); - - // compute electron and ion kinetic energy after emission - amrex::ParticleReal const u1sq_after = (u1x*u1x + u1y*u1y + u1z*u1z); - amrex::ParticleReal const u2sq_after = (u2x*u2x + u2y*u2y + u2z*u2z); - amrex::ParticleReal const gamma1_after = std::sqrt(1.0_prt + u1sq_after/c2); - amrex::ParticleReal const gamma2_after = std::sqrt(1.0_prt + u2sq_after/c2); - amrex::ParticleReal const KE1_after = weight1*m1*u1sq_after/(1.0_prt + gamma1_after); - amrex::ParticleReal const KE2_after = weight2*m2*u2sq_after/(1.0_prt + gamma2_after); - - // update energy lost to bremsstrahlung - amrex::ParticleReal const deltaE = KE1_before + KE2_before - KE1_after - KE2_after; - /* m_deltaE_bremsstrahlung += deltaE*weight1; // Joules */ - - if (m_create_photons) { - // compute photon energy in lab frame - amrex::ParticleReal const Ephoton_lab = deltaE/w_photon; // Joules - - p_mask[pair_index] = true; - p_pair_reaction_weight[pair_index] = w_photon; - p_product_data[pair_index] = Ephoton_lab; - } + // The electron and ion momentum are updated in the same place where the photon is created, + // in PhotonCreationFunc. + + p_product_data[pair_index] = Ephoton; + p_pair_reaction_weight[pair_index] = w_photon; + p_mask[pair_index] = 1; } @@ -461,7 +414,6 @@ public: bool m_computeSpeciesDensities = true; bool m_computeSpeciesTemperatures = false; - bool m_create_photons = true; bool m_need_product_data = true; amrex::ParticleReal m_multiplier; amrex::ParticleReal m_koT1_cut_default; diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp index 6ac5d4a4017..eec4b356627 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.cpp @@ -31,10 +31,6 @@ BremsstrahlungFunc::BremsstrahlungFunc (std::string const& collision_name, Multi WARPX_ALWAYS_ASSERT_WITH_MESSAGE(product_species.AmIA(), "BremsstrahlungFunc: The product species must be photons"); - bool create_photons = true; - pp_collision_name.query("create_photons", create_photons); - m_exe.m_create_photons = create_photons; - amrex::ParticleReal multiplier = 1._prt; pp_collision_name.query("multiplier", multiplier); WARPX_ALWAYS_ASSERT_WITH_MESSAGE(multiplier >= 1., diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H index 920460095e8..f75a7aab1bf 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H @@ -13,6 +13,7 @@ #include "Particles/ParticleCreation/SmartCopy.H" #include "Particles/MultiParticleContainer.H" #include "Particles/WarpXParticleContainer.H" +#include "Utils/ParticleUtils.H" #include #include @@ -90,17 +91,17 @@ public: AMREX_INLINE amrex::Vector operator() ( const index_type& n_total_pairs, - ParticleTileType& ptile1, ParticleTileType& /*ptile2*/, + ParticleTileType& ptile1, ParticleTileType& ptile2, const amrex::Vector& pc_products, ParticleTileType** AMREX_RESTRICT tile_products, - const amrex::ParticleReal& /*m1*/, const amrex::ParticleReal& /*m2*/, + const amrex::ParticleReal& m1, const amrex::ParticleReal& m2, const amrex::Vector& /*products_mass*/, const index_type* AMREX_RESTRICT p_mask, const amrex::Vector& products_np, const SmartCopy* AMREX_RESTRICT copy_species1, const SmartCopy* AMREX_RESTRICT /*copy_species2*/, const index_type* AMREX_RESTRICT p_pair_indices_1, - const index_type* AMREX_RESTRICT /*p_pair_indices_2*/, + const index_type* AMREX_RESTRICT p_pair_indices_2, const amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, const amrex::ParticleReal* AMREX_RESTRICT p_product_data ) const @@ -111,9 +112,15 @@ public: // Compute offset array and allocate memory for the produced species amrex::Gpu::DeviceVector offsets(n_total_pairs); - const auto total = amrex::Scan::ExclusiveSum(n_total_pairs, p_mask, offsets.data()); const index_type* AMREX_RESTRICT p_offsets = offsets.dataPtr(); + int total; + if (m_create_photons) { + total = amrex::Scan::ExclusiveSum(n_total_pairs, p_mask, offsets.data()); + } else { + total = 0; + } + amrex::Vector num_added_vec(m_num_product_species); for (int i = 0; i < m_num_product_species; i++) { @@ -123,6 +130,7 @@ public: } const auto soa_1 = ptile1.getParticleTileData(); + const auto soa_2 = ptile2.getParticleTileData(); // Create necessary GPU vectors, that will be used in the kernel below amrex::Vector soa_products; @@ -149,30 +157,92 @@ public: const index_type* AMREX_RESTRICT products_np_data = products_np.data(); #endif + bool create_photons = m_create_photons; + amrex::ParallelForRNG(n_total_pairs, [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept { if (p_mask[i]) { - const auto product_index = products_np_data[0] + p_offsets[i]; + constexpr auto c2 = PhysConst::c * PhysConst::c; + + int const i1 = static_cast(p_pair_indices_1[i]); + int const i2 = static_cast(p_pair_indices_2[i]); + + // Move to the ion rest frame and calculate the momentum and energy of the three particles + // This is done here to avoid having to save the photon momentum components + amrex::ParticleReal & u1x = soa_1.m_rdata[PIdx::ux][i1]; + amrex::ParticleReal & u1y = soa_1.m_rdata[PIdx::uy][i1]; + amrex::ParticleReal & u1z = soa_1.m_rdata[PIdx::uz][i1]; + amrex::ParticleReal & w1 = soa_1.m_rdata[PIdx::w][i1]; + amrex::ParticleReal & u2x = soa_2.m_rdata[PIdx::ux][i2]; + amrex::ParticleReal & u2y = soa_2.m_rdata[PIdx::uy][i2]; + amrex::ParticleReal & u2z = soa_2.m_rdata[PIdx::uz][i2]; + amrex::ParticleReal & w2 = soa_2.m_rdata[PIdx::w][i2]; + + amrex::ParticleReal u1x_rel = u1x; + amrex::ParticleReal u1y_rel = u1y; + amrex::ParticleReal u1z_rel = u1z; + ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, u2x, u2y, u2z); + + amrex::ParticleReal const u1sq_rel = (u1x_rel*u1x_rel + u1y_rel*u1y_rel + u1z_rel*u1z_rel); + amrex::ParticleReal const gamma1_rel = std::sqrt(1.0_prt + u1sq_rel/c2); - // Create product particle at position of particle 1 - copy_species1[0](soa_products_data[0], soa_1, static_cast(p_pair_indices_1[i]), - static_cast(product_index), engine); + amrex::ParticleReal const Ephoton = p_product_data[i]; + amrex::ParticleReal const wp = p_pair_reaction_weight[i]; - // Set the weight of the new particle to p_pair_reaction_weight[i] - soa_products_data[0].m_rdata[PIdx::w][product_index] = p_pair_reaction_weight[i]; + // Calculate electron and ion energy and momentum using conservation + amrex::ParticleReal const u1_rel = std::sqrt(u1sq_rel); + amrex::ParticleReal const A = u1_rel/PhysConst::c - Ephoton*wp/(w1*m1*c2); + amrex::ParticleReal const B = gamma1_rel + w2*m2/(w1*m1) - Ephoton*wp/(w1*m1*c2); + amrex::ParticleReal const D = w2*m2*w2*m2/(w1*m1*w1*m1) + A*A - B*B - 1.0_prt; + amrex::ParticleReal const a = 4.0_prt*(B*B - A*A); + amrex::ParticleReal const b = 4.0_prt*A*D; + amrex::ParticleReal const c = 4.0_prt*B*B - D*D; + // Use the "plus" solution since it keeps the electon motion forward + amrex::ParticleReal const u1_rel_after = (-b + std::sqrt(b*b - 4.0_prt*a*c))/(2.0_prt*a)*PhysConst::c; + amrex::ParticleReal const u2_rel_after = (A*PhysConst::c - u1_rel_after)*w1*m1/(w2*m2); - amrex::ParticleReal & ux1 = soa_products_data[0].m_rdata[PIdx::ux][product_index]; - amrex::ParticleReal & uy1 = soa_products_data[0].m_rdata[PIdx::uy][product_index]; - amrex::ParticleReal & uz1 = soa_products_data[0].m_rdata[PIdx::uz][product_index]; + u1x_rel *= u1_rel_after/u1_rel; + u1y_rel *= u1_rel_after/u1_rel; + u1z_rel *= u1_rel_after/u1_rel; + amrex::ParticleReal u2x_rel = u2_rel_after*u1x_rel/u1_rel_after; + amrex::ParticleReal u2y_rel = u2_rel_after*u1y_rel/u1_rel_after; + amrex::ParticleReal u2z_rel = u2_rel_after*u1z_rel/u1_rel_after; + amrex::ParticleReal p3x_rel = (Ephoton/PhysConst::c)*u1x_rel/u1_rel_after; + amrex::ParticleReal p3y_rel = (Ephoton/PhysConst::c)*u1y_rel/u1_rel_after; + amrex::ParticleReal p3z_rel = (Ephoton/PhysConst::c)*u1z_rel/u1_rel_after; - // Normalize out the electron velocity and multiply by the photon momentum, E/c - // Also, the photon momentum is normalized by m_e - amrex::ParticleReal const u1 = std::sqrt(ux1*ux1 + uy1*uy1 + uz1*uz1); - ux1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; - uy1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; - uz1 *= (p_product_data[i]/PhysConst::c)/u1/PhysConst::m_e; + // Lorentz transform electron and ion back to lab frame + ParticleUtils::doLorentzTransformWithU(u1x_rel, u1y_rel, u1z_rel, -u2x, -u2y, -u2z); + ParticleUtils::doLorentzTransformWithU(u2x_rel, u2y_rel, u2z_rel, -u2x, -u2y, -u2z); + ParticleUtils::doLorentzTransformWithP(p3x_rel, p3y_rel, p3z_rel, 0.0_prt, -u2x, -u2y, -u2z); + u1x = u1x_rel; + u1y = u1y_rel; + u1z = u1z_rel; + u2x = u2x_rel; + u2y = u2y_rel; + u2z = u2z_rel; + + if (create_photons) { + + // Create product particle at position of particle 1 + const auto product_index = products_np_data[0] + p_offsets[i]; + copy_species1[0](soa_products_data[0], soa_1, i1, + static_cast(product_index), engine); + + // Set the weight of the new particle to p_pair_reaction_weight[i] + soa_products_data[0].m_rdata[PIdx::w][product_index] = wp; + + amrex::ParticleReal & upx = soa_products_data[0].m_rdata[PIdx::ux][product_index]; + amrex::ParticleReal & upy = soa_products_data[0].m_rdata[PIdx::uy][product_index]; + amrex::ParticleReal & upz = soa_products_data[0].m_rdata[PIdx::uz][product_index]; + + // Lorentz transform to lab frame (from ion rest frame) + upx = p3x_rel/PhysConst::m_e; + upy = p3y_rel/PhysConst::m_e; + upz = p3z_rel/PhysConst::m_e; + } } }); @@ -206,6 +276,9 @@ private: // How many different type of species the collision produces int m_num_product_species; + // Whether photons are created + bool m_create_photons; + CollisionType m_collision_type; }; diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp index 120657a8382..46f1721a7e3 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.cpp @@ -32,4 +32,10 @@ PhotonCreationFunc::PhotonCreationFunc (const std::string& collision_name, WARPX_ABORT_WITH_MESSAGE("Unknown collision type in PhotonCreationFunc"); } + const amrex::ParmParse pp_collision_name(collision_name); + + bool create_photons = true; + pp_collision_name.query("create_photons", create_photons); + m_create_photons = create_photons; + } diff --git a/Source/Utils/ParticleUtils.H b/Source/Utils/ParticleUtils.H index 1423194bd72..fae255a4ad5 100644 --- a/Source/Utils/ParticleUtils.H +++ b/Source/Utils/ParticleUtils.H @@ -147,14 +147,52 @@ namespace ParticleUtils { amrex::ParticleReal const U = std::sqrt(Usq); if (U > 0._prt) { amrex::ParticleReal const gamma = std::sqrt(1._prt + Usq/c2); + // A nice way of calculating (gamma - 1) when it is small + amrex::ParticleReal const gammam1 = Usq/c2/(gamma + 1.0_prt); amrex::ParticleReal const Ux_hat = Ux/U; amrex::ParticleReal const Uy_hat = Uy/U; amrex::ParticleReal const Uz_hat = Uz/U; amrex::ParticleReal const P0 = std::sqrt(1._prt + (ux*ux + uy*uy + uz*uz)/c2); - amrex::ParticleReal const udotv = ux*Ux_hat + uy*Uy_hat + uz*Uz_hat; - ux = ux + (gamma - 1._prt)*udotv*Ux_hat - Ux*P0; - uy = uy + (gamma - 1._prt)*udotv*Uy_hat - Uy*P0; - uz = uz + (gamma - 1._prt)*udotv*Uz_hat - Uz*P0; + amrex::ParticleReal const udotn = ux*Ux_hat + uy*Uy_hat + uz*Uz_hat; + ux = ux + gammam1*udotn*Ux_hat - Ux*P0; + uy = uy + gammam1*udotn*Uy_hat - Uy*P0; + uz = uz + gammam1*udotn*Uz_hat - Uz*P0; + } + + } + + /** + * \brief Perform a Lorentz transformation of the given momentum + * to a frame moving with gamma*velocity (Ux, Uy, Uz) relative to the present one. + * + * @param[in,out] px,py,pz components of momentum vector in the current frame + * @param[in] Ux,Uy,Uz velocity of the new frame relative to the current one, + importantly these quantities are gamma * velocity + */ + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void doLorentzTransformWithP ( amrex::ParticleReal& px, amrex::ParticleReal& py, + amrex::ParticleReal& pz, amrex::ParticleReal mass, + amrex::ParticleReal const Ux, amrex::ParticleReal const Uy, + amrex::ParticleReal const Uz ) + { + using namespace amrex::literals; + + constexpr auto c2 = PhysConst::c * PhysConst::c; + + amrex::ParticleReal const Usq = Ux*Ux + Uy*Uy + Uz*Uz; + amrex::ParticleReal const U = std::sqrt(Usq); + if (U > 0._prt) { + amrex::ParticleReal const gamma = std::sqrt(1._prt + Usq/c2); + // A nice way of calculating (gamma - 1) when it is small + amrex::ParticleReal const gammam1 = Usq/c2/(gamma + 1.0_prt); + amrex::ParticleReal const Ux_hat = Ux/U; + amrex::ParticleReal const Uy_hat = Uy/U; + amrex::ParticleReal const Uz_hat = Uz/U; + amrex::ParticleReal const P0 = std::sqrt(px*px + py*py + pz*pz + mass*mass*c2); + amrex::ParticleReal const pdotn = px*Ux_hat + py*Uy_hat + pz*Uz_hat; + px = px + gammam1*pdotn*Ux_hat - Ux*P0/PhysConst::c; + py = py + gammam1*pdotn*Uy_hat - Uy*P0/PhysConst::c; + pz = pz + gammam1*pdotn*Uz_hat - Uz*P0/PhysConst::c; } } From 3e356fe9a5dcff919e0296dfa4669f9b4f151eaa Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 27 Jan 2025 09:32:12 -0800 Subject: [PATCH 36/38] Use const --- .../BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H index f75a7aab1bf..13632fe09fb 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H @@ -157,7 +157,7 @@ public: const index_type* AMREX_RESTRICT products_np_data = products_np.data(); #endif - bool create_photons = m_create_photons; + bool const create_photons = m_create_photons; amrex::ParallelForRNG(n_total_pairs, [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept @@ -174,11 +174,11 @@ public: amrex::ParticleReal & u1x = soa_1.m_rdata[PIdx::ux][i1]; amrex::ParticleReal & u1y = soa_1.m_rdata[PIdx::uy][i1]; amrex::ParticleReal & u1z = soa_1.m_rdata[PIdx::uz][i1]; - amrex::ParticleReal & w1 = soa_1.m_rdata[PIdx::w][i1]; + amrex::ParticleReal const & w1 = soa_1.m_rdata[PIdx::w][i1]; amrex::ParticleReal & u2x = soa_2.m_rdata[PIdx::ux][i2]; amrex::ParticleReal & u2y = soa_2.m_rdata[PIdx::uy][i2]; amrex::ParticleReal & u2z = soa_2.m_rdata[PIdx::uz][i2]; - amrex::ParticleReal & w2 = soa_2.m_rdata[PIdx::w][i2]; + amrex::ParticleReal const & w2 = soa_2.m_rdata[PIdx::w][i2]; amrex::ParticleReal u1x_rel = u1x; amrex::ParticleReal u1y_rel = u1y; From 0ed4e3ef26b705c3880e1414c4b3dd4fdf7b959e Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 27 Jan 2025 14:03:28 -0800 Subject: [PATCH 37/38] Expanding some expressions allows cancelling of some terms --- .../Bremsstrahlung/PhotonCreationFunc.H | 24 +++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H index 13632fe09fb..7c5e3a67960 100644 --- a/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/Bremsstrahlung/PhotonCreationFunc.H @@ -195,12 +195,32 @@ public: amrex::ParticleReal const u1_rel = std::sqrt(u1sq_rel); amrex::ParticleReal const A = u1_rel/PhysConst::c - Ephoton*wp/(w1*m1*c2); amrex::ParticleReal const B = gamma1_rel + w2*m2/(w1*m1) - Ephoton*wp/(w1*m1*c2); - amrex::ParticleReal const D = w2*m2*w2*m2/(w1*m1*w1*m1) + A*A - B*B - 1.0_prt; - amrex::ParticleReal const a = 4.0_prt*(B*B - A*A); + + // Expanding A**2 - B**2 cancels some terms and reduces the round off error + amrex::ParticleReal const AAmBB = -1.0_prt - 2.0_prt*gamma1_rel*w2*m2/(w1*m1) - w2*m2/(w1*m1)*w2*m2/(w1*m1) + +2.0_prt*(gamma1_rel + w2*m2/(w1*m1) - u1_rel/PhysConst::c)*Ephoton*wp/(w1*m1*c2); + amrex::ParticleReal const D = w2*m2*w2*m2/(w1*m1*w1*m1) + AAmBB - 1.0_prt; + + amrex::ParticleReal const a = 4.0_prt*(-AAmBB); amrex::ParticleReal const b = 4.0_prt*A*D; amrex::ParticleReal const c = 4.0_prt*B*B - D*D; + + // c could also be expanded but it doesn't seem to improve the round off + // u = u1_rel/PhysConst::c + // E = Ephoton*wp/(w1*m1*c2) + // w = w2*m2/(w1*m1) + // g = gamma1_rel + // cc1 = w*w*(- 4.0_prt*u*u + 8.0_prt*E*g - 4.0_prt*E*E) + // cc2 = w*(+ 8.0_prt*E*u*u - 8.0_prt*E*u*g + 8.0_prt*E + 8.0_prt*E*E*u - 8.0_prt*E*E*g) + // cc3 = (- 8.0_prt*E*u + 4.0_prt*u*u - 8.0_prt*E*E*u*u + 8.0_prt*E*E*u*g) + // c = cc1 + cc2 + cc3 + // Use the "plus" solution since it keeps the electon motion forward amrex::ParticleReal const u1_rel_after = (-b + std::sqrt(b*b - 4.0_prt*a*c))/(2.0_prt*a)*PhysConst::c; + + // The terms inside the sqrt can also be expanded but doesn't seem to improve the round off error + // u1_rel_after = (-b + 4.0_prt*std::sqrt(B*B*D*D + 4.0_prt*A*A*B*B - 4.0_prt*B*B*B*B))/(2.0_prt*a)*PhysConst::c + amrex::ParticleReal const u2_rel_after = (A*PhysConst::c - u1_rel_after)*w1*m1/(w2*m2); u1x_rel *= u1_rel_after/u1_rel; From 87107f20b0b059f2ce4d47028e9dc05febbb6462 Mon Sep 17 00:00:00 2001 From: David Grote Date: Mon, 27 Jan 2025 14:07:26 -0800 Subject: [PATCH 38/38] Update CI benchmark after adding conservation of momentum --- .../test_1d_collision_z_Bremsstrahlung.json | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json b/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json index 6adbbcc826f..0f006bba4e5 100644 --- a/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json +++ b/Regression/Checksum/benchmarks_json/test_1d_collision_z_Bremsstrahlung.json @@ -10,25 +10,25 @@ "jy": 0.0, "jz": 0.0 }, - "ions": { - "particle_momentum_x": 6.244058619001046e-18, - "particle_momentum_y": 6.266839027572808e-18, - "particle_momentum_z": 6.237935835816471e-18, - "particle_position_x": 0.0511999988707094, - "particle_weight": 5.470000000000001e+25 - }, "electrons": { - "particle_momentum_x": 1.8547285435878255e-26, - "particle_momentum_y": 1.8336044974219876e-26, - "particle_momentum_z": 7.780403913830049e-17, - "particle_position_x": 0.05119965192677569, + "particle_momentum_x": 1.7754689276894523e-26, + "particle_momentum_y": 1.8019937989491576e-26, + "particle_momentum_z": 7.780447460298363e-17, + "particle_position_x": 0.05119966127601259, "particle_weight": 5.470000000000001e+24 }, + "ions": { + "particle_momentum_x": 6.244058618052987e-18, + "particle_momentum_y": 6.266839026611909e-18, + "particle_momentum_z": 6.237937377390828e-18, + "particle_position_x": 0.051199998872927935, + "particle_weight": 5.470000000000001e+25 + }, "photons": { - "particle_momentum_x": 4.200876982473409e-27, - "particle_momentum_y": 4.172037413327287e-27, - "particle_momentum_z": 1.3125096507793955e-18, - "particle_position_x": 0.013402699474072983, - "particle_weight": 1.4387061523437505e+22 + "particle_momentum_x": 8.902000149229396e-25, + "particle_momentum_y": 9.02263331497841e-25, + "particle_momentum_z": 1.2715332613329798e-18, + "particle_position_x": 0.014344002048765016, + "particle_weight": 1.5380635742187505e+22 } }