Change alpha particle momentum initialization (#5033)

* Change alpha particle momentum initialization

Now alpha particles will be initializes with momentum calculated
via the two step alpha_1 channel, which has a higher cross section
than the alpha_0 channel currently used.

* Fix: Corrected Be mass

Corrected the mass of the Be to account for the mass defect due to being
an excited state. The conservation of energy tests now pass.

* Edit fusion and decay energy

* Modified released in fusion and decay steps of pB reaction. Ensures
the delta mc2 = energy released, at both the fusion and decay step

* Modified tolerence in pB fusion test

* Test one - abs. tol. in min energy of alpha2 and alpha3 is 20 keV
* Test two - abs. tol. in min energy of alpha2 and alpha3 is 200 keV

* Update checksum for tests

* Updated alpha particle energy checks

* Revert max/min energy check test

In `check_initial_energy2` the max and min sim. values are above and below theoretical calculations.
Tests now check variation is within an absolute tolerance.
Note: Total energy is still conserved.

* Increased atol for check_initial_energy2

* Changed max alpha particle energy check

Now the maximum alpha particle energy from simulation is compared to the
analytical value using relative tolerance, not abs tolerance.

* Fixed type

Author: PhysicsDan

Examples/Tests/nuclear_fusion/analysis_proton_boron_fusion.py

@@ -68,15 +68,15 @@
 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*scc.m_p # Beryllium 8 mass
+m_a = 4.00260325413*scc.m_u # Alpha (He4) mass
+m_be = (8.0053095729+0.00325283863)*scc.m_u # Be8* mass (3.03 MeV ex. state)
 Z_boron = 5.
 Z_proton = 1.
 E_Gamow = (Z_boron*Z_proton*np.pi*scc.fine_structure)**2*2.*m_reduced*scc.c**2
 E_Gamow_MeV = E_Gamow/MeV_to_Joule
 E_Gamow_keV = E_Gamow/keV_to_Joule
-E_fusion = 8.59009*MeV_to_Joule # Energy released during p + B -> alpha + Be
-E_decay = 0.0918984*MeV_to_Joule # Energy released during Be -> 2*alpha
+E_fusion = 5.55610759*MeV_to_Joule # Energy released during p + B -> alpha + Be*
+E_decay = 3.12600414*MeV_to_Joule # Energy released during Be* -> 2*alpha
 E_fusion_total = E_fusion + E_decay # Energy released during p + B -> 3*alpha
 
 ## Checks whether this is the 2D or the 3D test
@@ -493,10 +493,13 @@ def check_initial_energy1(data, E_com):
         energy_alpha3_simulation = energy_alpha_slice[4::6]
 
         assert(np.all(is_close(energy_alpha1_simulation, energy_alpha1_theory, rtol=5.e-8)))
-        assert(is_close(np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=1.e-2))
-        assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, rtol=1.e-2))
-        assert(is_close(np.amax(energy_alpha3_simulation), max_energy_alpha23, rtol=1.e-2))
-        assert(is_close(np.amin(energy_alpha3_simulation), min_energy_alpha23, rtol=1.e-2))
+        ## Check that the max / min value are comparable to the analytical value
+        ## The minimum value is checked to be within 20 keV of the analytical value
+        ## The maximum value is checked to be within 1% of the analytical value
+        assert(is_close(np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=1.e-2 ))
+        assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, atol=3.218e-15 ))
+        assert(is_close(np.amax(energy_alpha3_simulation), max_energy_alpha23, rtol=1.e-2 ))
+        assert(is_close(np.amin(energy_alpha3_simulation), min_energy_alpha23, atol=3.218e-15 ))
 
 def check_initial_energy2(data):
     ## In WarpX, the initial momentum of the alphas is computed assuming that the fusion process
@@ -578,13 +581,16 @@ def check_initial_energy2(data):
 
         assert(is_close(np.amax(energy_alpha1_simulation), max_energy_alpha1, rtol=1.e-2))
         assert(is_close(np.amin(energy_alpha1_simulation), min_energy_alpha1, rtol=1.e-2))
-        ## Tolerance is quite high below because we don't have a lot of alphas to produce good
+        ## Check that the max / min value are comparable to the analytical value
+        ## The minimum value is checked to be within 200 keV of the analytical value
+        ## The maximum value is checked to be within 5% of the analytical value
+        ## Tolerance is quite high because we don't have a lot of alphas to produce good
         ## statistics and an event like alpha1 emitted exactly in direction of proton & alpha2
         ## emitted exactly in direction opposite to Beryllium is somewhat rare.
-        assert(is_close(np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=2.5e-1))
-        assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, rtol=2.5e-1))
-        assert(is_close(np.amax(energy_alpha3_simulation), max_energy_alpha23, rtol=2.5e-1))
-        assert(is_close(np.amin(energy_alpha3_simulation), min_energy_alpha23, rtol=2.5e-1))
+        assert(is_close(np.amax(energy_alpha2_simulation), max_energy_alpha23, rtol=5e-2 ))
+        assert(is_close(np.amin(energy_alpha2_simulation), min_energy_alpha23, atol=3.218e-14 ))
+        assert(is_close(np.amax(energy_alpha3_simulation), max_energy_alpha23, rtol=5e-2 ))
+        assert(is_close(np.amin(energy_alpha3_simulation), min_energy_alpha23, atol=3.218e-14 ))
 
 def check_xy_isotropy(data):
     ## Checks that the alpha particles are emitted isotropically in x and y

Regression/Checksum/benchmarks_json/Proton_Boron_Fusion_2D.json

@@ -2,22 +2,6 @@
   "lev=0": {
     "rho": 0.0
   },
-  "alpha2": {
-    "particle_momentum_x": 4.057984510682467e-15,
-    "particle_momentum_y": 4.104188139725855e-15,
-    "particle_momentum_z": 4.17858000090827e-15,
-    "particle_position_x": 408793.7905852193,
-    "particle_position_y": 861780.8020495367,
-    "particle_weight": 5.078061191951185e+18
-  },
-  "alpha3": {
-    "particle_momentum_x": 5.017656304003558e-16,
-    "particle_momentum_y": 4.935595075276182e-16,
-    "particle_momentum_z": 4.867133212376827e-16,
-    "particle_position_x": 52678.192400911765,
-    "particle_position_y": 105483.59950020742,
-    "particle_weight": 1.5413633830148085e+27
-  },
   "boron1": {
     "particle_momentum_x": 0.0,
     "particle_momentum_y": 0.0,
@@ -26,6 +10,22 @@
     "particle_position_y": 81921136.14476715,
     "particle_weight": 128.00000000000261
   },
+  "boron3": {
+    "particle_momentum_x": 9.277692671587846e-15,
+    "particle_momentum_y": 9.268409636965691e-15,
+    "particle_momentum_z": 9.279446607709548e-15,
+    "particle_position_x": 4096178.1664224654,
+    "particle_position_y": 8192499.7060386725,
+    "particle_weight": 6.399486212205656e+30
+  },
+  "alpha5": {
+    "particle_momentum_x": 2.3827075626988002e-14,
+    "particle_momentum_y": 2.388035811027714e-14,
+    "particle_momentum_z": 2.4040151761604078e-14,
+    "particle_position_x": 2457556.8571638423,
+    "particle_position_y": 4914659.635379322,
+    "particle_weight": 3.839999999999998e-19
+  },
   "boron2": {
     "particle_momentum_x": 0.0,
     "particle_momentum_y": 0.0,
@@ -34,53 +34,37 @@
     "particle_position_y": 819270.9858143466,
     "particle_weight": 1.2799999998307316e+28
   },
-  "boron5": {
-    "particle_momentum_x": 0.0,
-    "particle_momentum_y": 0.0,
-    "particle_momentum_z": 0.0,
-    "particle_position_x": 409547.3312927569,
-    "particle_position_y": 819118.5558814355,
-    "particle_weight": 127.99999999999999
-  },
   "alpha1": {
-    "particle_momentum_x": 4.691909092431811e-15,
-    "particle_momentum_y": 4.6836958163275755e-15,
-    "particle_momentum_z": 4.657203546376977e-15,
+    "particle_momentum_x": 4.758865177822356e-15,
+    "particle_momentum_y": 4.7625215232714745e-15,
+    "particle_momentum_z": 4.734412976565431e-15,
     "particle_position_x": 463465.27131109464,
     "particle_position_y": 978207.6061359186,
     "particle_weight": 4.977661460251435e-28
   },
-  "proton2": {
-    "particle_momentum_x": 0.0,
-    "particle_momentum_y": 0.0,
-    "particle_momentum_z": 2.3723248133690294e-14,
-    "particle_position_x": 4095630.6981353555,
-    "particle_position_y": 8192073.551798361,
-    "particle_weight": 1.2885512788807322e+28
-  },
-  "proton1": {
+  "proton5": {
     "particle_momentum_x": 0.0,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 2.524872467113344e-13,
-    "particle_position_x": 40960140.72983792,
-    "particle_position_y": 81919772.69310114,
-    "particle_weight": 128.00000000000261
+    "particle_momentum_z": 1.7586062624930794e-15,
+    "particle_position_x": 409638.2877618571,
+    "particle_position_y": 819101.3225783394,
+    "particle_weight": 1.2800000000000004e+37
   },
-  "alpha5": {
-    "particle_momentum_x": 2.3343485859070736e-14,
-    "particle_momentum_y": 2.3451128753701046e-14,
-    "particle_momentum_z": 2.3579462789062662e-14,
-    "particle_position_x": 2457556.8571638423,
-    "particle_position_y": 4914659.635379322,
-    "particle_weight": 3.839999999999998e-19
+  "proton3": {
+    "particle_momentum_x": 1.6847290893972251e-15,
+    "particle_momentum_y": 1.6827074502304075e-15,
+    "particle_momentum_z": 1.6802489646490975e-15,
+    "particle_position_x": 2457270.6999197667,
+    "particle_position_y": 4914315.665267942,
+    "particle_weight": 1.279486212205663e+30
   },
-  "boron3": {
-    "particle_momentum_x": 9.277692671587846e-15,
-    "particle_momentum_y": 9.268409636965691e-15,
-    "particle_momentum_z": 9.279446607709548e-15,
-    "particle_position_x": 4096178.1664224654,
-    "particle_position_y": 8192499.7060386725,
-    "particle_weight": 6.399486212205656e+30
+  "alpha2": {
+    "particle_momentum_x": 4.140109871048478e-15,
+    "particle_momentum_y": 4.159837925039718e-15,
+    "particle_momentum_z": 4.251680389225861e-15,
+    "particle_position_x": 408793.7905852193,
+    "particle_position_y": 861780.8020495367,
+    "particle_weight": 5.078061191951185e+18
   },
   "proton4": {
     "particle_momentum_x": 0.0,
@@ -90,28 +74,44 @@
     "particle_position_y": 819198.7077771134,
     "particle_weight": 1.2800000000000004e+37
   },
-  "proton3": {
-    "particle_momentum_x": 1.6847290893972251e-15,
-    "particle_momentum_y": 1.6827074502304075e-15,
-    "particle_momentum_z": 1.6802489646490975e-15,
-    "particle_position_x": 2457270.6999197667,
-    "particle_position_y": 4914315.665267942,
-    "particle_weight": 1.279486212205663e+30
+  "proton1": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 2.524872467113344e-13,
+    "particle_position_x": 40960140.72983792,
+    "particle_position_y": 81919772.69310114,
+    "particle_weight": 128.00000000000261
+  },
+  "boron5": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 409547.3312927569,
+    "particle_position_y": 819118.5558814355,
+    "particle_weight": 127.99999999999999
   },
   "alpha4": {
-    "particle_momentum_x": 2.338084461204216e-14,
-    "particle_momentum_y": 2.3436156778849828e-14,
-    "particle_momentum_z": 2.35535708386288e-14,
+    "particle_momentum_x": 2.3854178157605226e-14,
+    "particle_momentum_y": 2.3882457581100904e-14,
+    "particle_momentum_z": 2.403451456378969e-14,
     "particle_position_x": 2457367.4582781536,
     "particle_position_y": 4915112.044373056,
     "particle_weight": 384.0000000000002
   },
-  "proton5": {
+  "proton2": {
     "particle_momentum_x": 0.0,
     "particle_momentum_y": 0.0,
-    "particle_momentum_z": 1.7586062624930794e-15,
-    "particle_position_x": 409638.2877618571,
-    "particle_position_y": 819101.3225783394,
-    "particle_weight": 1.2800000000000004e+37
+    "particle_momentum_z": 2.3723248133690294e-14,
+    "particle_position_x": 4095630.6981353555,
+    "particle_position_y": 8192073.551798361,
+    "particle_weight": 1.2885512788807322e+28
+  },
+  "alpha3": {
+    "particle_momentum_x": 5.079136896829917e-16,
+    "particle_momentum_y": 5.075922501147803e-16,
+    "particle_momentum_z": 4.962151258214859e-16,
+    "particle_position_x": 52678.192400911765,
+    "particle_position_y": 105483.59950020742,
+    "particle_weight": 1.5413633830148085e+27
   }
 }