diff options
Diffstat (limited to 'Source/Particles')
7 files changed, 249 insertions, 32 deletions
diff --git a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp index fc70b39ca..7e78d95ed 100644 --- a/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp +++ b/Source/Particles/Collision/BinaryCollision/BinaryCollisionUtils.cpp @@ -26,9 +26,9 @@ namespace BinaryCollisionUtils{ auto& species1 = mypc->GetParticleContainerFromName(species_names[0]); auto& species2 = mypc->GetParticleContainerFromName(species_names[1]); - if ((species1.AmIA<PhysicalSpecies::hydrogen>() && species2.AmIA<PhysicalSpecies::boron>()) + if ((species1.AmIA<PhysicalSpecies::hydrogen>() && species2.AmIA<PhysicalSpecies::boron11>()) || - (species1.AmIA<PhysicalSpecies::boron>() && species2.AmIA<PhysicalSpecies::hydrogen>()) + (species1.AmIA<PhysicalSpecies::boron11>() && species2.AmIA<PhysicalSpecies::hydrogen>()) ) { return NuclearFusionType::ProtonBoron; diff --git a/Source/Particles/Collision/BinaryCollision/NuclearFusionFunc.H b/Source/Particles/Collision/BinaryCollision/NuclearFusionFunc.H index bf3c9118b..895254ed3 100644 --- a/Source/Particles/Collision/BinaryCollision/NuclearFusionFunc.H +++ b/Source/Particles/Collision/BinaryCollision/NuclearFusionFunc.H @@ -57,6 +57,10 @@ public: { 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); diff --git a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H index 15bfb44e5..33c4ec520 100644 --- a/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H +++ b/Source/Particles/Collision/BinaryCollision/ParticleCreationFunc.H @@ -214,7 +214,7 @@ public: p_num_products_device[0]; ProtonBoronFusionInitializeMomentum(soa_1, soa_2, soa_products_data[0], p_pair_indices_1[i], p_pair_indices_2[i], - product_start_index, m1, m2); + product_start_index, m1, m2, engine); } } }); diff --git a/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionCrossSection.H b/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionCrossSection.H index 330cceaf7..fad4d87de 100644 --- a/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionCrossSection.H +++ b/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionCrossSection.H @@ -105,5 +105,4 @@ amrex::ParticleReal ProtonBoronFusionCrossSection (const amrex::ParticleReal& E_ return cross_section_b*barn_to_sqm; } - #endif // PROTON_BORON_FUSION_CROSS_SECTION_H diff --git a/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionInitializeMomentum.H b/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionInitializeMomentum.H index 511c2d637..415555193 100644 --- a/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionInitializeMomentum.H +++ b/Source/Particles/Collision/BinaryCollision/ProtonBoronFusionInitializeMomentum.H @@ -9,10 +9,16 @@ #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; @@ -22,35 +28,227 @@ namespace { /** * \brief This function initializes the momentum of the alpha particles produced from - * proton-boron fusion. It needs to be implemented and currently only zero-initializes the - * momenta. + * 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*/, - const SoaData_type soa_alpha, - const index_type& /*idx_1*/, const index_type& /*idx_2*/, + 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::ParticleReal& m1, const amrex::ParticleReal& m2, + const amrex::RandomEngine& engine) { - soa_alpha.m_rdata[PIdx::ux][idx_alpha_start] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uy][idx_alpha_start] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uz][idx_alpha_start] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 1] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 1] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 1] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 2] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 2] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 2] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 3] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 3] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 3] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 4] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 4] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 4] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 5] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 5] = amrex::ParticleReal(0.); - soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 5] = amrex::ParticleReal(0.); + // 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; } } diff --git a/Source/Particles/Collision/BinaryCollision/SingleNuclearFusionEvent.H b/Source/Particles/Collision/BinaryCollision/SingleNuclearFusionEvent.H index 28030467d..a1f3db324 100644 --- a/Source/Particles/Collision/BinaryCollision/SingleNuclearFusionEvent.H +++ b/Source/Particles/Collision/BinaryCollision/SingleNuclearFusionEvent.H @@ -92,9 +92,9 @@ void SingleNuclearFusionEvent (const amrex::ParticleReal& u1x, const amrex::Part 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); + 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; diff --git a/Source/Particles/SpeciesPhysicalProperties.H b/Source/Particles/SpeciesPhysicalProperties.H index 6a3fd12b6..b70f95c50 100644 --- a/Source/Particles/SpeciesPhysicalProperties.H +++ b/Source/Particles/SpeciesPhysicalProperties.H @@ -18,7 +18,7 @@ #include <string> enum struct PhysicalSpecies{unspecified=0, electron, positron, photon, hydrogen, helium, boron, - carbon, nitrogen, oxygen, argon, copper, xenon}; + boron10, boron11, carbon, nitrogen, oxygen, argon, copper, xenon}; namespace species { @@ -43,6 +43,10 @@ namespace species return PhysicalSpecies::helium; if( species=="boron" ) return PhysicalSpecies::boron; + if( species=="boron10" ) + return PhysicalSpecies::boron10; + if( species=="boron11" ) + return PhysicalSpecies::boron11; if( species=="carbon" ) return PhysicalSpecies::carbon; if( species=="nitrogen" ) @@ -77,6 +81,10 @@ namespace species return PhysConst::q_e * amrex::Real(2.0); case PhysicalSpecies::boron: return PhysConst::q_e * amrex::Real(5.0); + case PhysicalSpecies::boron10: + return PhysConst::q_e * amrex::Real(5.0); + case PhysicalSpecies::boron11: + return PhysConst::q_e * amrex::Real(5.0); case PhysicalSpecies::carbon: return PhysConst::q_e * amrex::Real(6.0); case PhysicalSpecies::nitrogen: @@ -113,6 +121,10 @@ namespace species return PhysConst::m_p * amrex::Real(3.97369); case PhysicalSpecies::boron: return PhysConst::m_p * amrex::Real(10.7319); + case PhysicalSpecies::boron10: + return PhysConst::m_p * amrex::Real(9.94060); + case PhysicalSpecies::boron11: + return PhysConst::m_p * amrex::Real(10.9298); case PhysicalSpecies::carbon: return PhysConst::m_e * amrex::Real(22032.0); case PhysicalSpecies::nitrogen: @@ -149,6 +161,10 @@ namespace species return "helium"; case PhysicalSpecies::boron: return "boron"; + case PhysicalSpecies::boron10: + return "boron10"; + case PhysicalSpecies::boron11: + return "boron11"; case PhysicalSpecies::carbon: return "carbon"; case PhysicalSpecies::nitrogen: |