Use more consistent species types in fusion module (ECP-WarpX#4480)

* Force user to select helium4 in the p-B fusion reaction

* Update helium4 mass in fusion module and tests

* Replace proton with hydrogen1

* Update hydrogen1 mass in fusion module and tests

* Update benchmarks

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.002602*scc.m_u # Alpha mass
-m_be = 7.94748*m_p # Beryllium 8 mass
+m_a = 4.00260325413*scc.m_u # Alpha 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,11 +33,12 @@ 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
 
-proton1.species_type = proton
+proton1.species_type = hydrogen1
 proton1.injection_style = "NRandomPerCell"
 proton1.num_particles_per_cell = 80000
 proton1.profile = constant
@@ -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
 
@@ -63,14 +64,14 @@ boron1.momentum_function_uz(x,y,z) = -sqrt(2*m_reduced*Energy_step*(floor(z)**2)
 boron1.do_not_push = 1
 boron1.do_not_deposit = 1
 
-alpha1.species_type = helium
+alpha1.species_type = helium4
 alpha1.do_not_push = 1
 alpha1.do_not_deposit = 1
 
 my_constants.background_dens = 1.e26
 my_constants.beam_dens = 1.e20
 
-proton2.species_type = proton
+proton2.species_type = hydrogen1
 proton2.injection_style = "NRandomPerCell"
 proton2.num_particles_per_cell = 8000
 proton2.profile = "parse_density_function"
@@ -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
 
@@ -94,19 +95,19 @@ boron2.momentum_distribution_type = "constant"
 boron2.do_not_push = 1
 boron2.do_not_deposit = 1
 
-alpha2.species_type = helium
+alpha2.species_type = helium4
 alpha2.do_not_push = 1
 alpha2.do_not_deposit = 1
 
 my_constants.temperature = 44. * keV_to_J
 
-proton3.species_type = proton
+proton3.species_type = hydrogen1
 proton3.injection_style = "NRandomPerCell"
 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
 
@@ -120,19 +121,19 @@ boron3.theta = temperature/(m_b11*clight**2)
 boron3.do_not_push = 1
 boron3.do_not_deposit = 1
 
-alpha3.species_type = helium
+alpha3.species_type = helium4
 alpha3.do_not_push = 1
 alpha3.do_not_deposit = 1
 
 my_constants.proton4_energy = 550*keV_to_J
 
-proton4.species_type = proton
+proton4.species_type = hydrogen1
 proton4.injection_style = "NRandomPerCell"
 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
 
@@ -145,17 +146,17 @@ boron4.momentum_distribution_type = "constant"
 boron4.do_not_push = 1
 boron4.do_not_deposit = 1
 
-alpha4.species_type = helium
+alpha4.species_type = helium4
 alpha4.do_not_push = 1
 alpha4.do_not_deposit = 1
 
-proton5.species_type = proton
+proton5.species_type = hydrogen1
 proton5.injection_style = "NRandomPerCell"
 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
 
@@ -168,7 +169,7 @@ boron5.momentum_distribution_type = "constant"
 boron5.do_not_push = 1
 boron5.do_not_deposit = 1
 
-alpha5.species_type = helium
+alpha5.species_type = helium4
 alpha5.do_not_push = 1
 alpha5.do_not_deposit = 1
 

Examples/Tests/nuclear_fusion/inputs_proton_boron_3d

@@ -33,11 +33,12 @@ 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
 
-proton1.species_type = proton
+proton1.species_type = hydrogen1
 proton1.injection_style = "NRandomPerCell"
 proton1.num_particles_per_cell = 10000
 proton1.profile = constant
@@ -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
 
@@ -63,14 +64,14 @@ boron1.momentum_function_uz(x,y,z) = -sqrt(2*m_reduced*Energy_step*(floor(z)**2)
 boron1.do_not_push = 1
 boron1.do_not_deposit = 1
 
-alpha1.species_type = helium
+alpha1.species_type = helium4
 alpha1.do_not_push = 1
 alpha1.do_not_deposit = 1
 
 my_constants.background_dens = 1.e26
 my_constants.beam_dens = 1.e20
 
-proton2.species_type = proton
+proton2.species_type = hydrogen1
 proton2.injection_style = "NRandomPerCell"
 proton2.num_particles_per_cell = 1000
 proton2.profile = "parse_density_function"
@@ -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
 
@@ -94,19 +95,19 @@ boron2.momentum_distribution_type = "constant"
 boron2.do_not_push = 1
 boron2.do_not_deposit = 1
 
-alpha2.species_type = helium
+alpha2.species_type = helium4
 alpha2.do_not_push = 1
 alpha2.do_not_deposit = 1
 
 my_constants.temperature = 44. * keV_to_J
 
-proton3.species_type = proton
+proton3.species_type = hydrogen1
 proton3.injection_style = "NRandomPerCell"
 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
 
@@ -120,19 +121,19 @@ boron3.theta = temperature/(m_b11*clight**2)
 boron3.do_not_push = 1
 boron3.do_not_deposit = 1
 
-alpha3.species_type = helium
+alpha3.species_type = helium4
 alpha3.do_not_push = 1
 alpha3.do_not_deposit = 1
 
 my_constants.proton4_energy = 550*keV_to_J
 
-proton4.species_type = proton
+proton4.species_type = hydrogen1
 proton4.injection_style = "NRandomPerCell"
 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
 
@@ -145,17 +146,17 @@ boron4.momentum_distribution_type = "constant"
 boron4.do_not_push = 1
 boron4.do_not_deposit = 1
 
-alpha4.species_type = helium
+alpha4.species_type = helium4
 alpha4.do_not_push = 1
 alpha4.do_not_deposit = 1
 
-proton5.species_type = proton
+proton5.species_type = hydrogen1
 proton5.injection_style = "NRandomPerCell"
 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
 
@@ -168,7 +169,7 @@ boron5.momentum_distribution_type = "constant"
 boron5.do_not_push = 1
 boron5.do_not_deposit = 1
 
-alpha5.species_type = helium
+alpha5.species_type = helium4
 alpha5.do_not_push = 1
 alpha5.do_not_deposit = 1