Update hydrogen1 mass in fusion module and tests

Examples/Tests/nuclear_fusion/analysis_proton_boron_fusion.py

@@ -65,11 +65,11 @@
 keV_to_Joule = scc.e*1e3
 MeV_to_Joule = scc.e*1e6
 barn_to_square_meter = 1.e-28
-m_p = scc.m_p # Proton mass
+m_p = 1.00782503223*scc.m_u # Proton mass
 m_b = 11.00930536*scc.m_u # Boron 11 mass
 m_reduced = m_p*m_b/(m_p+m_b)
 m_a = 4.00260325413*scc.m_u # Alpha mass
-m_be = 7.94748*m_p # Beryllium 8 mass
+m_be = 7.94748*scc.m_p # Beryllium 8 mass
 Z_boron = 5.
 Z_proton = 1.
 E_Gamow = (Z_boron*Z_proton*np.pi*scc.fine_structure)**2*2.*m_reduced*scc.c**2
@@ -395,7 +395,7 @@ def p_sq_boron_frame_to_E_COM_frame(p_proton_sq):
     # Use invariant E**2 - p**2c**2 of 4-momentum norm to compute energy in center of mass frame
     E_com = np.sqrt(E_lab**2 - p_proton_sq*scc.c**2)
     # Corresponding kinetic energy
-    E_com_kin = E_com - (m_b+scc.m_p)*scc.c**2
+    E_com_kin = E_com - (m_b+m_p)*scc.c**2
     return E_com_kin*(p_proton_sq>0.)
 
 def p_sq_to_kinetic_energy(p_sq, m):
@@ -525,7 +525,7 @@ def check_initial_energy2(data):
         ## Proton kinetic energy in the lab frame before fusion
         E_proton_nonrelativistic = Energy_step*slice_number**2
         ## Corresponding square norm of proton momentum
-        p_proton_sq = 2.*scc.m_p*E_proton_nonrelativistic
+        p_proton_sq = 2.*m_p*E_proton_nonrelativistic
         ## Kinetic energy in the lab frame after
         ## proton + boron 11 -> alpha + beryllium 8
         E_after_fusion = E_proton_nonrelativistic + E_fusion

Examples/Tests/nuclear_fusion/inputs_proton_boron_2d

@@ -33,7 +33,8 @@ particles.species_names = proton1 boron1 alpha1 proton2 boron2 alpha2 proton3 bo
                           proton4 boron4 alpha4 proton5 boron5 alpha5
 
 my_constants.m_b11 =  11.00930536*m_u # Boron 11 mass
-my_constants.m_reduced = m_p*m_b11/(m_p+m_b11)
+my_constants.m_h = 1.00782503223*m_u # Hydrogen 1 mass
+my_constants.m_reduced = m_h*m_b11/(m_h+m_b11)
 my_constants.keV_to_J = 1.e3*q_e
 my_constants.Energy_step = 22. * keV_to_J
 
@@ -46,7 +47,7 @@ proton1.momentum_distribution_type = "parse_momentum_function"
 proton1.momentum_function_ux(x,y,z) = 0.
 proton1.momentum_function_uy(x,y,z) = 0.
 ## Thanks to the floor, all particles in the same cell have the exact same momentum
-proton1.momentum_function_uz(x,y,z) = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_p*clight)
+proton1.momentum_function_uz(x,y,z) = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_h*clight)
 proton1.do_not_push = 1
 proton1.do_not_deposit = 1
 
@@ -81,7 +82,7 @@ proton2.momentum_distribution_type = "parse_momentum_function"
 proton2.momentum_function_ux(x,y,z) = 0.
 proton2.momentum_function_uy(x,y,z) = 0.
 proton2.momentum_function_uz(x,y,z) = "if(x - floor(x) < 0.1,
-                                          0., sqrt(2*m_p*Energy_step*(floor(z)**2))/(m_p*clight))"
+                                          0., sqrt(2*m_h*Energy_step*(floor(z)**2))/(m_h*clight))"
 proton2.do_not_push = 1
 proton2.do_not_deposit = 1
 
@@ -106,7 +107,7 @@ proton3.num_particles_per_cell = 4800
 proton3.profile = constant
 proton3.density = 1.e28
 proton3.momentum_distribution_type = "maxwell_boltzmann"
-proton3.theta = temperature/(m_p*clight**2)
+proton3.theta = temperature/(m_h*clight**2)
 proton3.do_not_push = 1
 proton3.do_not_deposit = 1
 
@@ -132,7 +133,7 @@ proton4.num_particles_per_cell = 800
 proton4.profile = "constant"
 proton4.density = 1.e35
 proton4.momentum_distribution_type = "constant"
-proton4.uz = sqrt(2*m_p*proton4_energy)/(m_p*clight)
+proton4.uz = sqrt(2*m_h*proton4_energy)/(m_h*clight)
 proton4.do_not_push = 1
 proton4.do_not_deposit = 1
 
@@ -155,7 +156,7 @@ proton5.num_particles_per_cell = 800
 proton5.profile = "constant"
 proton5.density = 1.e35
 proton5.momentum_distribution_type = "constant"
-proton5.uz = sqrt(2*m_p*proton4_energy)/(m_p*clight)
+proton5.uz = sqrt(2*m_h*proton4_energy)/(m_h*clight)
 proton5.do_not_push = 1
 proton5.do_not_deposit = 1
 

Examples/Tests/nuclear_fusion/inputs_proton_boron_3d

@@ -33,7 +33,8 @@ particles.species_names = proton1 boron1 alpha1 proton2 boron2 alpha2 proton3 bo
                           proton4 boron4 alpha4 proton5 boron5 alpha5
 
 my_constants.m_b11 = 11.00930536*m_u # Boron 11 mass
-my_constants.m_reduced = m_p*m_b11/(m_p+m_b11)
+my_constants.m_h = 1.00782503223*m_u # Hydrogen 1 mass
+my_constants.m_reduced = m_h*m_b11/(m_h+m_b11)
 my_constants.keV_to_J = 1.e3*q_e
 my_constants.Energy_step = 22. * keV_to_J
 
@@ -46,7 +47,7 @@ proton1.momentum_distribution_type = "parse_momentum_function"
 proton1.momentum_function_ux(x,y,z) = 0.
 proton1.momentum_function_uy(x,y,z) = 0.
 ## Thanks to the floor, all particles in the same cell have the exact same momentum
-proton1.momentum_function_uz(x,y,z) = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_p*clight)
+proton1.momentum_function_uz(x,y,z) = sqrt(2*m_reduced*Energy_step*(floor(z)**2))/(m_h*clight)
 proton1.do_not_push = 1
 proton1.do_not_deposit = 1
 
@@ -81,7 +82,7 @@ proton2.momentum_distribution_type = "parse_momentum_function"
 proton2.momentum_function_ux(x,y,z) = 0.
 proton2.momentum_function_uy(x,y,z) = 0.
 proton2.momentum_function_uz(x,y,z) = "if(y - floor(y) < 0.1,
-                                          0., sqrt(2*m_p*Energy_step*(floor(z)**2))/(m_p*clight))"
+                                          0., sqrt(2*m_h*Energy_step*(floor(z)**2))/(m_h*clight))"
 proton2.do_not_push = 1
 proton2.do_not_deposit = 1
 
@@ -106,7 +107,7 @@ proton3.num_particles_per_cell = 600
 proton3.profile = constant
 proton3.density = 1.e28
 proton3.momentum_distribution_type = "maxwell_boltzmann"
-proton3.theta = temperature/(m_p*clight**2)
+proton3.theta = temperature/(m_h*clight**2)
 proton3.do_not_push = 1
 proton3.do_not_deposit = 1
 
@@ -132,7 +133,7 @@ proton4.num_particles_per_cell = 100
 proton4.profile = "constant"
 proton4.density = 1.e35
 proton4.momentum_distribution_type = "constant"
-proton4.uz = sqrt(2*m_p*proton4_energy)/(m_p*clight)
+proton4.uz = sqrt(2*m_h*proton4_energy)/(m_h*clight)
 proton4.do_not_push = 1
 proton4.do_not_deposit = 1
 
@@ -155,7 +156,7 @@ proton5.num_particles_per_cell = 100
 proton5.profile = "constant"
 proton5.density = 1.e35
 proton5.momentum_distribution_type = "constant"
-proton5.uz = sqrt(2*m_p*proton4_energy)/(m_p*clight)
+proton5.uz = sqrt(2*m_h*proton4_energy)/(m_h*clight)
 proton5.do_not_push = 1
 proton5.do_not_deposit = 1
 

Source/Particles/Collision/BinaryCollision/NuclearFusion/ProtonBoronFusionCrossSection.H

@@ -42,8 +42,9 @@ amrex::ParticleReal ProtonBoronFusionCrossSectionNevins (const amrex::ParticleRe
     // 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 m_boron = 11.00930536_prt * PhysConst::m_u;
+    constexpr amrex::ParticleReal m_hydrogen = 1.00782503223 * PhysConst::m_u;
+    constexpr amrex::ParticleReal m_reduced = m_boron / (one_pr + m_boron/m_hydrogen);
     constexpr amrex::ParticleReal gamow_factor = m_reduced / 2._prt *
                         (PhysConst::q_e*PhysConst::q_e * Z_boron /
                            (2._prt*PhysConst::ep0*PhysConst::hbar)) *