diff options
Diffstat (limited to 'Source/Particles/Collision/BinaryCollision/NuclearFusion')
4 files changed, 824 insertions, 0 deletions
diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H new file mode 100644 index 000000000..47575f42c --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H @@ -0,0 +1,224 @@ +/* Copyright 2021 Neil Zaim + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef NUCLEAR_FUSION_FUNC_H_ +#define NUCLEAR_FUSION_FUNC_H_ + +#include "SingleNuclearFusionEvent.H" + +#include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H" +#include "Particles/Pusher/GetAndSetPosition.H" +#include "Particles/MultiParticleContainer.H" +#include "Particles/WarpXParticleContainer.H" +#include "Utils/WarpXUtil.H" +#include "WarpX.H" + +#include <AMReX_Algorithm.H> +#include <AMReX_DenseBins.H> +#include <AMReX_ParmParse.H> +#include <AMReX_Random.H> +#include <AMReX_REAL.H> +#include <AMReX_Vector.H> + +/** + * \brief This functor does binary nuclear fusions on a single cell. + * Particles of the two reacting species are paired with each other and for each pair we compute + * if a fusion event occurs. If so, we fill a mask (input parameter p_mask) with true so that + * product particles corresponding to a given pair can be effectively created in the particle + * creation functor. + * This functor also reads and contains the fusion multiplier. + */ +class NuclearFusionFunc{ + // Define shortcuts for frequently-used type names + using ParticleType = WarpXParticleContainer::ParticleType; + using ParticleBins = amrex::DenseBins<ParticleType>; + using index_type = ParticleBins::index_type; + using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + +public: + /** + * \brief Default constructor of the NuclearFusionFunc class. + */ + NuclearFusionFunc () = default; + + /** + * \brief Constructor of the NuclearFusionFunc 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 + */ + NuclearFusionFunc (const std::string collision_name, MultiParticleContainer const * const mypc, + const bool isSameSpecies) : m_isSameSpecies(isSameSpecies) + { + using namespace amrex::literals; + +#ifdef AMREX_SINGLE_PRECISION_PARTICLES + amrex::Abort("Nuclear fusion module does not currently work with single precision"); +#endif + + m_fusion_type = BinaryCollisionUtils::get_nuclear_fusion_type(collision_name, mypc); + + amrex::ParmParse pp_collision_name(collision_name); + amrex::Vector<std::string> product_species_name; + pp_collision_name.getarr("product_species", product_species_name); + + if (m_fusion_type == NuclearFusionType::ProtonBoron) + { + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + product_species_name.size() == 1, + "ERROR: Proton-boron must contain exactly one product species"); + auto& product_species = mypc->GetParticleContainerFromName(product_species_name[0]); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + product_species.AmIA<PhysicalSpecies::helium>(), + "ERROR: Product species of proton-boron fusion must be of type helium"); + } + + // default fusion multiplier + m_fusion_multiplier = 1.0_rt; + queryWithParser(pp_collision_name, "fusion_multiplier", m_fusion_multiplier); + // default fusion probability threshold + m_probability_threshold = 0.02_rt; + queryWithParser(pp_collision_name, "fusion_probability_threshold", + m_probability_threshold); + // default fusion probability target_value + m_probability_target_value = 0.002_rt; + queryWithParser(pp_collision_name, "fusion_probability_target_value", + m_probability_target_value); + } + + /** + * \brief operator() of the NuclearFusionFunc class. Performs nuclear fusions at the cell level + * using the algorithm described in Higginson et al., Journal of Computational Physics 388, + * 439-453 (2019). Note that this function does not yet create the product particles, but + * instead fills an array p_mask that stores which collisions result in a fusion event. + * + * Also note that there are three main differences between this implementation and the + * algorithm described in Higginson's paper: + * - 1) The transformation from the lab frame to the center of mass frame is nonrelativistic + * in Higginson's paper. Here, we implement a relativistic generalization. + * - 2) The behaviour when the estimated fusion probability is greater than one is not + * specified in Higginson's paper. Here, we provide an implementation using two runtime + * dependent parameters (fusion probability threshold and fusion probability target value). See + * documentation for more details. + * - 3) Here, we divide the weight of a particle by the number of times it is paired with other + * particles. This was not explicitly specified in Higginson's paper. + * + * @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] 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 if a fusion event occurs for a given + * pair. It is only needed here to store information that will be used later on when actually + * creating the product particles. + * @param[out] p_pair_indices_1,p_pair_indices_2 arrays that store the indices of the + * particles of a given pair. They are only needed here to store information that will be used + * later on when actually creating the product particles. + * @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[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* I1, index_type* I2, + SoaData_type soa_1, SoaData_type soa_2, + GetParticlePosition /*get_position_1*/, GetParticlePosition /*get_position_2*/, + 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 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::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 int NI1 = I1e - I1s; + const int NI2 = I2e - I2s; + const int max_N = amrex::max(NI1,NI2); + + int i1 = I1s; + int i2 = I2s; + int pair_index = cell_start_pair; + + // Because the number of particles of each species is not always equal (NI1 != NI2 + // in general), some macroparticles will be paired with multiple macroparticles of the + // other species and we need to decrease their weight accordingly. + // c1 corresponds to the minimum number of times a particle of species 1 will be paired + // with a particle of species 2. Same for c2. + const int c1 = amrex::max(NI2/NI1,1); + const int c2 = amrex::max(NI1/NI2,1); + + // multiplier ratio to take into account unsampled pairs + int multiplier_ratio; + if (m_isSameSpecies) + { + multiplier_ratio = 2*max_N - 1; + } + else + { + multiplier_ratio = max_N; + } + + for (int k = 0; k < max_N; ++k) + { + // c1k : how many times the current particle of species 1 is paired with a particle + // of species 2. Same for c2k. + const int c1k = (k%NI1 < max_N%NI1) ? c1 + 1: c1; + const int c2k = (k%NI2 < max_N%NI2) ? c2 + 1: c2; + SingleNuclearFusionEvent( + u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], + u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], + m1, m2, w1[ I1[i1] ]/c1k, w2[ I2[i2] ]/c2k, + dt, dV, pair_index, p_mask, p_pair_reaction_weight, + m_fusion_multiplier, multiplier_ratio, + m_probability_threshold, + m_probability_target_value, + m_fusion_type, engine); + p_pair_indices_1[pair_index] = I1[i1]; + p_pair_indices_2[pair_index] = I2[i2]; + ++i1; if ( i1 == static_cast<int>(I1e) ) { i1 = I1s; } + ++i2; if ( i2 == static_cast<int>(I2e) ) { i2 = I2s; } + ++pair_index; + } + + } + +private: + // Factor used to increase the number of fusion reaction by decreasing the weight of the + // produced particles + amrex::Real m_fusion_multiplier; + // If the fusion multiplier is too high and results in a fusion probability that approaches + // 1, there is a risk of underestimating the total fusion yield. In these cases, we reduce + // the fusion multiplier used in a given collision. m_probability_threshold is the fusion + // probability threshold above which we reduce the fusion multiplier. + // m_probability_target_value is the target probability used to determine by how much + // the fusion multiplier should be reduced. + amrex::Real m_probability_threshold; + amrex::Real m_probability_target_value; + NuclearFusionType m_fusion_type; + bool m_isSameSpecies; +}; + +#endif // NUCLEAR_FUSION_FUNC_H_ diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H new file mode 100644 index 000000000..ee81a620a --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H @@ -0,0 +1,154 @@ +/* Copyright 2021 Neil Zaim + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef PROTON_BORON_FUSION_CROSS_SECTION_H +#define PROTON_BORON_FUSION_CROSS_SECTION_H + +#include "Utils/WarpXConst.H" + +#include <AMReX_REAL.H> + +#include <cmath> + +/** + * \brief Computes the total proton-boron fusion cross section in the range 0 < E < 3.5 MeV using + * the analytical fits given in W.M. Nevins and R. Swain, Nuclear Fusion, 40, 865 (2000). + * For the record, note that there is a typo in equation (1) of this paper: the total cross section + * should read S(E)/E*exp(-sqrt(E_G/E)) instead of S(E)/E*exp(sqrt(E_G/E)) (minus sign in the + * exponential). + * + * @param[in] E_keV the kinetic energy of the proton-boron pair in its center of mass frame, in + * keV. + * @return The total cross section in barn. + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::ParticleReal ProtonBoronFusionCrossSectionNevins (const amrex::ParticleReal& E_keV) +{ + using namespace amrex::literals; + + // If kinetic energy is 0, return a 0 cross section and avoid later division by 0. + if (E_keV == 0._prt) {return 0._prt;} + + // Fits also use energy in MeV + const amrex::ParticleReal E_MeV = E_keV*1.e-3_prt; + constexpr amrex::ParticleReal joule_to_MeV = 1.e-6_prt/PhysConst::q_e; + + // Compute Gamow factor, in MeV + constexpr auto one_pr = 1._prt; + constexpr auto Z_boron = 5._prt; + constexpr amrex::ParticleReal m_boron = 10.7319_prt * PhysConst::m_p; + constexpr amrex::ParticleReal m_reduced = m_boron / (one_pr + m_boron/PhysConst::m_p); + constexpr amrex::ParticleReal gamow_factor = m_reduced / 2._prt * + (PhysConst::q_e*PhysConst::q_e * Z_boron / + (2._prt*PhysConst::ep0*PhysConst::hbar)) * + (PhysConst::q_e*PhysConst::q_e * Z_boron / + (2._prt*PhysConst::ep0*PhysConst::hbar)) * + joule_to_MeV; + + // Compute astrophysical factor, in MeV barn, using the fits + constexpr auto E_lim1 = 400._prt; // Limits between the different fit regions + constexpr auto E_lim2 = 642._prt; + amrex::ParticleReal astrophysical_factor; + if (E_keV < E_lim1) + { + constexpr auto C0 = 197._prt; + constexpr auto C1 = 0.24_prt; + constexpr auto C2 = 2.31e-4_prt; + constexpr auto AL = 1.82e4_prt; + constexpr auto EL = 148._prt; + constexpr auto dEL_sq = 2.35_prt*2.35_prt; + astrophysical_factor = C0 + C1*E_keV + C2*E_keV*E_keV + + AL/((E_keV - EL)*(E_keV - EL) + dEL_sq); + } + else if (E_keV < E_lim2) + { + constexpr auto D0 = 330._prt; + constexpr auto D1 = 66.1_prt; + constexpr auto D2 = -20.3_prt; + constexpr auto D5 = -1.58_prt; + const amrex::ParticleReal E_norm = (E_keV-400._prt) * 1.e-2_prt; + astrophysical_factor = D0 + D1*E_norm + D2*E_norm*E_norm + D5*std::pow(E_norm,5); + } + else + { + constexpr auto A0 = 2.57e6_prt; + constexpr auto A1 = 5.67e5_prt; + constexpr auto A2 = 1.34e5_prt; + constexpr auto A3 = 5.68e5_prt; + constexpr auto E0 = 581.3_prt; + constexpr auto E1 = 1083._prt; + constexpr auto E2 = 2405._prt; + constexpr auto E3 = 3344._prt; + constexpr auto dE0_sq = 85.7_prt*85.7_prt; + constexpr auto dE1_sq = 234._prt*234._prt; + constexpr auto dE2_sq = 138._prt*138._prt; + constexpr auto dE3_sq = 309._prt*309._prt; + constexpr auto B = 4.38_prt; + astrophysical_factor = A0 / ((E_keV-E0)*(E_keV-E0) + dE0_sq) + + A1 / ((E_keV-E1)*(E_keV-E1) + dE1_sq) + + A2 / ((E_keV-E2)*(E_keV-E2) + dE2_sq) + + A3 / ((E_keV-E3)*(E_keV-E3) + dE3_sq) + B; + } + + // Compute cross section, in barn + return astrophysical_factor/E_MeV*std::exp(-std::sqrt(gamow_factor/E_MeV)); +} + +/** + * \brief Computes the total proton-boron fusion cross section in the range E > 3.5 MeV using a + * simple power law fit of the data presented in Buck et al., Nuclear Physics A, 398(2), 189-202 + * (1983) (data can also be found in the EXFOR database). + * + * @param[in] E_keV the kinetic energy of the proton-boron pair in its center of mass frame, in + * keV. + * @return The total cross section in barn. + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::ParticleReal ProtonBoronFusionCrossSectionBuck (const amrex::ParticleReal& E_keV) +{ + using namespace amrex::literals; + + constexpr amrex::ParticleReal E_start_fit = 3500._prt; // Fit starts at 3.5 MeV + // cross section at E = E_start_fit, in barn + constexpr amrex::ParticleReal cross_section_start_fit = 0.2168440845211521_prt; + constexpr amrex::ParticleReal slope_fit = -2.661840717596765; + + // Compute fitted value + return cross_section_start_fit*std::pow(E_keV/E_start_fit, slope_fit); +} + +/** + * \brief Computes the total proton-boron fusion cross section. When E_kin_star < 3.5 MeV, we use + * the analytical fits given in W.M. Nevins and R. Swain, Nuclear Fusion, 40, 865 (2000). When + * E_kin_star > 3.5 MeV, we use a simple power law fit of the data presented in Buck et al., + * Nuclear Physics A, 398(2), 189-202 (1983). Both fits return the same value for + * E_kin_star = 3.5 MeV. + * + * @param[in] E_kin_star the kinetic energy of the proton-boron pair in its center of mass frame, + * in SI units. + * @return The total cross section in SI units (square meters). + */ +AMREX_GPU_HOST_DEVICE AMREX_INLINE +amrex::ParticleReal ProtonBoronFusionCrossSection (const amrex::ParticleReal& E_kin_star) +{ + using namespace amrex::literals; + + // Fits use energy in keV + constexpr amrex::ParticleReal joule_to_keV = 1.e-3_prt/PhysConst::q_e; + const amrex::ParticleReal E_keV = E_kin_star*joule_to_keV; + constexpr amrex::ParticleReal E_threshold = 3500._prt; + + const amrex::ParticleReal cross_section_b = (E_keV <= E_threshold) ? + ProtonBoronFusionCrossSectionNevins(E_keV) : + ProtonBoronFusionCrossSectionBuck(E_keV); + + // Convert cross section to SI units: barn to square meter + constexpr auto barn_to_sqm = 1.e-28_prt; + return cross_section_b*barn_to_sqm; +} + +#endif // PROTON_BORON_FUSION_CROSS_SECTION_H diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H new file mode 100644 index 000000000..415555193 --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionInitializeMomentum.H @@ -0,0 +1,256 @@ +/* Copyright 2021 Neil Zaim + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H +#define PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H + +#include "Particles/WarpXParticleContainer.H" +#include "Utils/ParticleUtils.H" +#include "Utils/WarpXConst.H" + +#include <AMReX_DenseBins.H> +#include <AMReX_Random.H> +#include <AMReX_REAL.H> + +#include <cmath> +#include <limits> + +namespace { + // Define shortcuts for frequently-used type names + using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType; + using ParticleType = WarpXParticleContainer::ParticleType; + using ParticleBins = amrex::DenseBins<ParticleType>; + using index_type = ParticleBins::index_type; + + /** + * \brief This function initializes the momentum of the alpha particles produced from + * proton-boron fusion. The momentum is initialized by assuming that the fusion of a proton + * with a boron nucleus into 3 alphas takes place in two steps. In the first step, the proton + * and the boron fuse into a beryllium nucleus and an alpha particle. In the second step, the + * beryllium decays into two alpha particles. The first step produces 8.59009 MeV of kinetic + * energy while the second step produces 91.8984 keV of kinetic energy. This two-step process + * is considered to be the dominant process of proton+boron fusion into alphas (see + * Becker et al., Zeitschrift für Physik A Atomic Nuclei, 327(3), 341-355 (1987)). + * For each step, we assume in this function that the particles are emitted isotropically in + * the corresponding center of mass frame (center of mass frame of proton + boron for the + * creation of first alpha+beryllium and rest frame of beryllium for the creation of second and + * third alphas). This isotropic assumption is exact for the second step but is only an + * approximation for the first step. + * + * @param[in] soa_1 struct of array data of the first colliding species (can be either proton + * or boron) + * @param[in] soa_2 struct of array data of the second colliding species (can be either proton + * or boron) + * @param[out] soa_alpha struct of array data of the alpha species + * @param[in] idx_1 index of first colliding macroparticle + * @param[in] idx_2 index of second colliding macroparticle + * @param[in] idx_alpha_start index of first produced alpha macroparticle + * @param[in] m1 mass of first colliding species + * @param[in] m2 mass of second colliding species + * @param[in] engine the random engine + */ + AMREX_GPU_HOST_DEVICE AMREX_INLINE + void ProtonBoronFusionInitializeMomentum ( + const SoaData_type& soa_1, const SoaData_type& soa_2, + SoaData_type& soa_alpha, + const index_type& idx_1, const index_type& idx_2, + const index_type& idx_alpha_start, + const amrex::ParticleReal& m1, const amrex::ParticleReal& m2, + const amrex::RandomEngine& engine) + { + // General notations in this function: + // x_sq denotes the square of x + // x_star denotes the value of x in the proton+boron center of mass frame + // x_Bestar denotes the value of x in the Beryllium rest frame + + using namespace amrex::literals; + + constexpr amrex::ParticleReal mev_to_joule = PhysConst::q_e*1.e6_prt; + // Energy produced in the fusion reaction proton + boron11 -> Beryllium8 + alpha + // cf. Janis book of proton-induced cross-sections (2019) + constexpr amrex::ParticleReal E_fusion = 8.59009_prt*mev_to_joule; + // Energy produced when Beryllium8 decays into two alphas + // cf. JEFF-3.3 radioactive decay data library (2017) + constexpr amrex::ParticleReal E_decay = 0.0918984_prt*mev_to_joule; + + constexpr amrex::ParticleReal m_alpha = PhysConst::m_p * 3.97369_prt; + constexpr amrex::ParticleReal m_beryllium = PhysConst::m_p * 7.94748_prt; + constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c; + constexpr amrex::ParticleReal inv_csq = 1._prt / ( c_sq ); + // Rest energy of proton+boron + const amrex::ParticleReal E_rest_pb = (m1 + m2)*c_sq; + // Rest energy of alpha+beryllium + constexpr amrex::ParticleReal E_rest_abe = (m_alpha + m_beryllium)*c_sq; + + // These constexprs underflow in single precision because we use SI units and sometimes the + // code won't compile. Nuclear fusion module does not currently work with single precision + // anyways so these #ifdef are just here to make sure that the code always compiles with + // single precision. +#ifdef AMREX_SINGLE_PRECISION_PARTICLES + const amrex::ParticleReal ma_sq = m_alpha*m_alpha; + const amrex::ParticleReal mBe_sq = m_beryllium*m_beryllium; +#else + constexpr amrex::ParticleReal ma_sq = m_alpha*m_alpha; + constexpr amrex::ParticleReal mBe_sq = m_beryllium*m_beryllium; +#endif + + // Normalized momentum of colliding particles (proton and boron) + const amrex::ParticleReal u1x = soa_1.m_rdata[PIdx::ux][idx_1]; + const amrex::ParticleReal u1y = soa_1.m_rdata[PIdx::uy][idx_1]; + const amrex::ParticleReal u1z = soa_1.m_rdata[PIdx::uz][idx_1]; + const amrex::ParticleReal u2x = soa_2.m_rdata[PIdx::ux][idx_2]; + const amrex::ParticleReal u2y = soa_2.m_rdata[PIdx::uy][idx_2]; + const amrex::ParticleReal u2z = soa_2.m_rdata[PIdx::uz][idx_2]; + + // Compute Lorentz factor gamma in the lab frame + const amrex::ParticleReal g1 = std::sqrt( 1._prt + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_csq ); + const amrex::ParticleReal g2 = std::sqrt( 1._prt + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_csq ); + + // Compute momenta + const amrex::ParticleReal p1x = u1x * m1; + const amrex::ParticleReal p1y = u1y * m1; + const amrex::ParticleReal p1z = u1z * m1; + const amrex::ParticleReal p2x = u2x * m2; + const amrex::ParticleReal p2y = u2y * m2; + const amrex::ParticleReal p2z = u2z * m2; + // Square norm of the total (sum between the two particles) momenta in the lab frame + const amrex::ParticleReal p_total_sq = std::pow(p1x+p2x, 2) + + std::pow(p1y+p2y, 2) + + std::pow(p1z+p2z, 2); + + // Total energy of proton+boron in the lab frame + const amrex::ParticleReal E_lab = (m1 * g1 + m2 * g2) * c_sq; + // Total energy squared of proton+boron in the center of mass frame, calculated using the + // Lorentz invariance of the four-momentum norm + const amrex::ParticleReal E_star_sq = E_lab*E_lab - c_sq*p_total_sq; + // Total energy squared of beryllium+alpha in the center of mass frame + // In principle, the term - E_rest_pb + E_rest_abe + E_fusion is not needed and equal to + // zero (i.e. the energy liberated during fusion is equal to the mass difference). However, + // due to possible inconsistencies in how the mass is defined in the code (e.g. currently, + // the mass of hydrogen is the mass of the proton, not including the electron, while the + // mass of the other elements is the atomic mass, that includes the electrons mass), it is + // probably more robust to subtract the rest masses and to add the fusion energy to the + // total kinetic energy. + const amrex::ParticleReal E_star_f_sq = std::pow(std::sqrt(E_star_sq) + - E_rest_pb + E_rest_abe + E_fusion, 2); + + // Square of the norm of the momentum of Beryllium or alpha in the center of mass frame + // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle + const amrex::ParticleReal p_star_f_sq = + E_star_f_sq*0.25_prt*inv_csq - (ma_sq + mBe_sq)*c_sq*0.5_prt + + std::pow(c_sq, 3)*0.25_prt*std::pow(ma_sq - mBe_sq, 2)/E_star_f_sq; + + // Compute momentum of first alpha in the center of mass frame, assuming isotropic + // distribution + amrex::ParticleReal px_star, py_star, pz_star; + ParticleUtils::RandomizeVelocity(px_star, py_star, pz_star, std::sqrt(p_star_f_sq), + engine); + + // Next step is to convert momentum of first alpha to lab frame + amrex::ParticleReal px_alpha1, py_alpha1, pz_alpha1; + // Preliminary calculation: compute center of mass velocity vc + const amrex::ParticleReal mass_g = m1 * g1 + m2 * g2; + const amrex::ParticleReal vcx = (p1x+p2x) / mass_g; + const amrex::ParticleReal vcy = (p1y+p2y) / mass_g; + const amrex::ParticleReal vcz = (p1z+p2z) / mass_g; + const amrex::ParticleReal vc_sq = vcx*vcx + vcy*vcy + vcz*vcz; + + // Convert momentum of first alpha to lab frame, using equation (13) of F. Perez et al., + // Phys.Plasmas.19.083104 (2012) + if ( vc_sq > std::numeric_limits<amrex::ParticleReal>::min() ) + { + const amrex::ParticleReal gc = 1._prt / std::sqrt( 1._prt - vc_sq*inv_csq ); + const amrex::ParticleReal g_star = std::sqrt(1._prt + p_star_f_sq / (ma_sq*c_sq)); + const amrex::ParticleReal vcDps = vcx*px_star + vcy*py_star + vcz*pz_star; + const amrex::ParticleReal factor0 = (gc-1._prt)/vc_sq; + const amrex::ParticleReal factor = factor0*vcDps + m_alpha*g_star*gc; + px_alpha1 = px_star + vcx * factor; + py_alpha1 = py_star + vcy * factor; + pz_alpha1 = pz_star + vcz * factor; + } + else // If center of mass velocity is zero, we are already in the lab frame + { + px_alpha1 = px_star; + py_alpha1 = py_star; + pz_alpha1 = pz_star; + } + + // Compute momentum of beryllium in lab frame, using total momentum conservation + const amrex::ParticleReal px_Be = p1x + p2x - px_alpha1; + const amrex::ParticleReal py_Be = p1y + p2y - py_alpha1; + const amrex::ParticleReal pz_Be = p1z + p2z - pz_alpha1; + + // Compute momentum norm of second and third alphas in Beryllium rest frame + // Factor 0.5 is here because each alpha only gets half of the decay energy + constexpr amrex::ParticleReal gamma_Bestar = (1._prt + 0.5_prt*E_decay/(m_alpha*c_sq)); + constexpr amrex::ParticleReal gamma_Bestar_sq_minus_one = gamma_Bestar*gamma_Bestar - 1._prt; + const amrex::ParticleReal p_Bestar = m_alpha*PhysConst::c*std::sqrt(gamma_Bestar_sq_minus_one); + + // Compute momentum of second alpha in Beryllium rest frame, assuming isotropic distribution + amrex::ParticleReal px_Bestar, py_Bestar, pz_Bestar; + ParticleUtils::RandomizeVelocity(px_Bestar, py_Bestar, pz_Bestar, p_Bestar, engine); + + // Next step is to convert momentum of second alpha to lab frame + amrex::ParticleReal px_alpha2, py_alpha2, pz_alpha2; + // Preliminary calculation: compute Beryllium velocity v_Be + const amrex::ParticleReal p_Be_sq = px_Be*px_Be + py_Be*py_Be + pz_Be*pz_Be; + const amrex::ParticleReal g_Be = std::sqrt(1._prt + p_Be_sq / (mBe_sq*c_sq)); + const amrex::ParticleReal mg_Be = m_beryllium*g_Be; + const amrex::ParticleReal v_Bex = px_Be / mg_Be; + const amrex::ParticleReal v_Bey = py_Be / mg_Be; + const amrex::ParticleReal v_Bez = pz_Be / mg_Be; + const amrex::ParticleReal v_Be_sq = v_Bex*v_Bex + v_Bey*v_Bey + v_Bez*v_Bez; + + // Convert momentum of second alpha to lab frame, using equation (13) of F. Perez et al., + // Phys.Plasmas.19.083104 (2012) + if ( v_Be_sq > std::numeric_limits<amrex::ParticleReal>::min() ) + { + const amrex::ParticleReal vcDps = v_Bex*px_Bestar + v_Bey*py_Bestar + v_Bez*pz_Bestar; + const amrex::ParticleReal factor0 = (g_Be-1._prt)/v_Be_sq; + const amrex::ParticleReal factor = factor0*vcDps + m_alpha*gamma_Bestar*g_Be; + px_alpha2 = px_Bestar + v_Bex * factor; + py_alpha2 = py_Bestar + v_Bey * factor; + pz_alpha2 = pz_Bestar + v_Bez * factor; + } + else // If Beryllium velocity is zero, we are already in the lab frame + { + px_alpha2 = px_Bestar; + py_alpha2 = py_Bestar; + pz_alpha2 = pz_Bestar; + } + + // Compute momentum of third alpha in lab frame, using total momentum conservation + const amrex::ParticleReal px_alpha3 = px_Be - px_alpha2; + const amrex::ParticleReal py_alpha3 = py_Be - py_alpha2; + const amrex::ParticleReal pz_alpha3 = pz_Be - pz_alpha2; + + // Fill alpha species momentum data with the computed momentum (note that we actually + // create 6 alphas, 3 at the position of the proton and 3 at the position of the boron, so + // each computed momentum is used twice) + soa_alpha.m_rdata[PIdx::ux][idx_alpha_start] = px_alpha1/m_alpha; + soa_alpha.m_rdata[PIdx::uy][idx_alpha_start] = py_alpha1/m_alpha; + soa_alpha.m_rdata[PIdx::uz][idx_alpha_start] = pz_alpha1/m_alpha; + soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 1] = px_alpha1/m_alpha; + soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 1] = py_alpha1/m_alpha; + soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 1] = pz_alpha1/m_alpha; + soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 2] = px_alpha2/m_alpha; + soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 2] = py_alpha2/m_alpha; + soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 2] = pz_alpha2/m_alpha; + soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 3] = px_alpha2/m_alpha; + soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 3] = py_alpha2/m_alpha; + soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 3] = pz_alpha2/m_alpha; + soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 4] = px_alpha3/m_alpha; + soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 4] = py_alpha3/m_alpha; + soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 4] = pz_alpha3/m_alpha; + soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 5] = px_alpha3/m_alpha; + soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 5] = py_alpha3/m_alpha; + soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 5] = pz_alpha3/m_alpha; + } + +} + +#endif // PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H new file mode 100644 index 000000000..c8633b52b --- /dev/null +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusion/SingleNuclearFusionEvent.H @@ -0,0 +1,190 @@ +/* Copyright 2021 Neil Zaim + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + */ + +#ifndef SINGLE_NUCLEAR_FUSION_EVENT_H_ +#define SINGLE_NUCLEAR_FUSION_EVENT_H_ + +#include "ProtonBoronFusionCrossSection.H" + +#include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H" +#include "Utils/WarpXConst.H" + +#include <AMReX_Algorithm.H> +#include <AMReX_Random.H> +#include <AMReX_REAL.H> + +#include <cmath> + + +/** + * \brief This function computes whether the collision between two particles result in a + * nuclear fusion event, using the algorithm described in Higginson et al., Journal of + * Computational Physics 388, 439-453 (2019). If nuclear fusion occurs, the mask is set to true + * for that given pair of particles and the weight of the produced particles is stored in + * p_pair_reaction_weight. + * + * @tparam index_type type of the index argument + * @param[in] u1x,u1y,u1z momenta of the first colliding particle + * @param[in] u2x,u2y,u2z momenta of the second colliding particle + * @param[in] m1,m2 masses + * @param[in] w1,w2 effective weight of the colliding particles + * @param[in] dt is the time step length between two collision calls. + * @param[in] dV is the volume of the corresponding cell. + * @param[in] pair_index is the index of the colliding pair + * @param[out] p_mask is a mask that will be set to true if fusion occurs for that pair + * @param[out] p_pair_reaction_weight stores the weight of the product particles + * @param[in] fusion_multiplier factor used to increase the number of fusion events by + * decreasing the weight of the produced particles + * @param[in] multiplier_ratio factor used to take into account unsampled pairs (i.e. the fact + * that a particle only collides with one or few particles of the other species) + * @param[in] probability_threshold probability threshold above which we decrease the fusion + * multiplier + * @param[in] probability_target_value if the probability threshold is exceeded, this is used + * to determine by how much the fusion multiplier is reduced + * @param[in] fusion_type the physical fusion process to model + * @param[in] engine the random engine. + */ +template <typename index_type> +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::ParticleReal& u1y, + const amrex::ParticleReal& u1z, const amrex::ParticleReal& u2x, + const amrex::ParticleReal& u2y, const amrex::ParticleReal& u2z, + const amrex::ParticleReal& m1, const amrex::ParticleReal& m2, + amrex::ParticleReal w1, amrex::ParticleReal w2, + const amrex::Real& dt, const amrex::Real& dV, const int& pair_index, + index_type* AMREX_RESTRICT p_mask, + amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight, + const amrex::Real& fusion_multiplier, + const int& multiplier_ratio, + const amrex::Real& probability_threshold, + const amrex::Real& probability_target_value, + const NuclearFusionType& fusion_type, + const amrex::RandomEngine& engine) +{ + // General notations in this function: + // x_sq denotes the square of x + // x_star denotes the value of x in the center of mass frame + + const amrex::ParticleReal w_min = amrex::min(w1, w2); + const amrex::ParticleReal w_max = amrex::max(w1, w2); + + constexpr auto one_pr = amrex::ParticleReal(1.); + constexpr auto inv_two_pr = amrex::ParticleReal(1./2.); + constexpr auto inv_four_pr = amrex::ParticleReal(1./4.); + constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c; + constexpr amrex::ParticleReal inv_csq = one_pr / ( c_sq ); + + const amrex::ParticleReal m1_sq = m1*m1; + const amrex::ParticleReal m2_sq = m2*m2; + + // Compute Lorentz factor gamma in the lab frame + const amrex::ParticleReal g1 = std::sqrt( one_pr + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_csq ); + const amrex::ParticleReal g2 = std::sqrt( one_pr + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_csq ); + + // Compute momenta + const amrex::ParticleReal p1x = u1x * m1; + const amrex::ParticleReal p1y = u1y * m1; + const amrex::ParticleReal p1z = u1z * m1; + const amrex::ParticleReal p2x = u2x * m2; + const amrex::ParticleReal p2y = u2y * m2; + const amrex::ParticleReal p2z = u2z * m2; + // Square norm of the total (sum between the two particles) momenta in the lab frame + const amrex::ParticleReal p_total_sq = std::pow(p1x+p2x, 2) + + std::pow(p1y+p2y, 2) + + std::pow(p1z+p2z, 2); + + // Total energy in the lab frame + const amrex::ParticleReal E_lab = (m1 * g1 + m2 * g2) * c_sq; + // Total energy squared in the center of mass frame, calculated using the Lorentz invariance + // of the four-momentum norm + const amrex::ParticleReal E_star_sq = E_lab*E_lab - c_sq*p_total_sq; + + // Kinetic energy in the center of mass frame + const amrex::ParticleReal E_kin_star = std::sqrt(E_star_sq) - (m1 + m2)*c_sq; + + // Compute fusion cross section as a function of kinetic energy in the center of mass frame + auto fusion_cross_section = amrex::ParticleReal(0.); + if (fusion_type == NuclearFusionType::ProtonBoron) + { + fusion_cross_section = ProtonBoronFusionCrossSection(E_kin_star); + } + + // Square of the norm of the momentum of one of the particles in the center of mass frame + // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle + const amrex::ParticleReal p_star_sq = + E_star_sq*inv_four_pr*inv_csq - (m1_sq + m2_sq)*c_sq*inv_two_pr + + std::pow(c_sq, 3)*inv_four_pr*std::pow(m1_sq - m2_sq, 2)/E_star_sq; + + // Lorentz factors in the center of mass frame + const amrex::ParticleReal g1_star = std::sqrt(1 + p_star_sq / (m1_sq*c_sq)); + const amrex::ParticleReal g2_star = std::sqrt(1 + p_star_sq / (m2_sq*c_sq)); + + // relative velocity in the center of mass frame + const amrex::ParticleReal v_rel = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) + + one_pr/(m2*g2_star)); + + // Fusion cross section and relative velocity are computed in the center of mass frame. + // On the other hand, the particle densities (weight over volume) in the lab frame are used. To + // take into account this discrepancy, we need to multiply the fusion probability by the ratio + // between the Lorentz factors in the COM frame and the Lorentz factors in the lab frame + // (see Perez et al., Phys.Plasmas.19.083104 (2012)) + const amrex::ParticleReal lab_to_COM_factor = g1_star*g2_star/(g1*g2); + + // First estimate of probability to have fusion reaction + amrex::ParticleReal probability_estimate = multiplier_ratio * fusion_multiplier * + lab_to_COM_factor * w_max * fusion_cross_section * v_rel * dt / dV; + + // Effective fusion multiplier + amrex::ParticleReal fusion_multiplier_eff = fusion_multiplier; + + // If the fusion probability is too high and the fusion multiplier greater than one, we risk to + // systematically underestimate the fusion yield. In this case, we reduce the fusion multiplier + // to reduce the fusion probability + if (probability_estimate > probability_threshold) + { + // We aim for a fusion probability of probability_target_value but take into account + // the constraint that the fusion_multiplier cannot be smaller than one + fusion_multiplier_eff = amrex::max(fusion_multiplier * + probability_target_value / probability_estimate , one_pr); + probability_estimate *= fusion_multiplier_eff/fusion_multiplier; + } + + // Compute actual fusion probability that is always between zero and one + // In principle this is obtained by computing 1 - exp(-probability_estimate) + // However, the computation of this quantity can fail numerically when probability_estimate is + // too small (e.g. exp(-probability_estimate) returns 1 and the computation returns 0). + // In this case, we simply use "probability_estimate" instead of 1 - exp(-probability_estimate) + // The threshold exp_threshold at which we switch between the two formulas is determined by the + // fact that computing the exponential is only useful if it can resolve the x^2/2 term of its + // Taylor expansion, i.e. the square of probability_estimate should be greater than the + // machine epsilon. +#ifdef AMREX_SINGLE_PRECISION_PARTICLES + constexpr auto exp_threshold = amrex::ParticleReal(1.e-3); +#else + constexpr auto exp_threshold = amrex::ParticleReal(5.e-8); +#endif + const amrex::ParticleReal probability = (probability_estimate < exp_threshold) ? + probability_estimate: one_pr - std::exp(-probability_estimate); + + // Get a random number + amrex::ParticleReal random_number = amrex::Random(engine); + + // If we have a fusion event, set the mask the true and fill the product weight array + if (random_number < probability) + { + p_mask[pair_index] = true; + p_pair_reaction_weight[pair_index] = w_min/fusion_multiplier_eff; + } + else + { + p_mask[pair_index] = false; + } + +} + + +#endif // SINGLE_NUCLEAR_FUSION_EVENT_H_ |