Prevent underflow in calculation of speed (#2811)

openmc-dev · Jan 3, 2024 · b0926ea · b0926ea
1 parent 2e4811b
commit b0926ea
Show file tree

Hide file tree

Showing 9 changed files with 69 additions and 60 deletions.
diff --git a/src/particle.cpp b/src/particle.cpp
@@ -59,11 +59,20 @@ double Particle::speed() const
     break;
   }
 
-  // Calculate inverse of Lorentz factor
-  const double inv_gamma = mass / (this->E() + mass);
+  if (this->E() < 1.0e-9 * mass) {
+    // If the energy is much smaller than the mass, revert to non-relativistic
+    // formula. The 1e-9 criterion is specifically chosen as the point below
+    // which the error from using the non-relativistic formula is less than the
+    // round-off eror when using the relativistic formula (see analysis at
+    // https://gist.github.com/paulromano/da3b473fe3df33de94b265bdff0c7817)
+    return C_LIGHT * std::sqrt(2 * this->E() / mass);
+  } else {
+    // Calculate inverse of Lorentz factor
+    const double inv_gamma = mass / (this->E() + mass);
 
-  // Calculate speed via v = c * sqrt(1 - γ^-2)
-  return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma);
+    // Calculate speed via v = c * sqrt(1 - γ^-2)
+    return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma);
+  }
 }
 
 void Particle::move_distance(double length)

diff --git a/tests/regression_tests/mgxs_library_condense/results_true.dat b/tests/regression_tests/mgxs_library_condense/results_true.dat
@@ -168,10 +168,10 @@
 3      2  2  1         1   total  1.0  0.130701
   mesh 1       group in nuclide          mean     std. dev.
        x  y  z                                             
-0      1  1  1        1   total  5.304284e-07  2.560239e-08
-2      1  2  1        1   total  4.940320e-07  2.410416e-08
-1      2  1  1        1   total  5.587366e-07  3.382787e-08
-3      2  2  1        1   total  5.232209e-07  2.482800e-08
+0      1  1  1        1   total  5.304285e-07  2.560239e-08
+2      1  2  1        1   total  4.940321e-07  2.410417e-08
+1      2  1  1        1   total  5.587365e-07  3.382787e-08
+3      2  2  1        1   total  5.232210e-07  2.482800e-08
   mesh 1       group in nuclide      mean std. dev.
        x  y  z                                     
 0      1  1  1        1   total  0.026032  0.001172

diff --git a/tests/regression_tests/mgxs_library_distribcell/results_true.dat b/tests/regression_tests/mgxs_library_distribcell/results_true.dat
@@ -51,7 +51,7 @@
                                                                                                                                                                                                                                                                                                                                                                                                  sum(distribcell)  group out nuclide  mean  std. dev.
 0  ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),)          1   total   1.0   0.084331
                                                                                                                                                                                                                                                                                                                                                                                                  sum(distribcell)  group in nuclide          mean     std. dev.
-0  ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),)         1   total  5.253873e-07  2.168461e-08
+0  ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),)         1   total  5.253873e-07  2.168462e-08
                                                                                                                                                                                                                                                                                                                                                                                                  sum(distribcell)  group in nuclide      mean  std. dev.
 0  ((0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, ...),)         1   total  0.094147   0.003701
                                                                                                                                                                                                                                                                                                                                                                                                  sum(distribcell)  group in  group out nuclide      mean  std. dev.

diff --git a/tests/regression_tests/mgxs_library_hdf5/results_true.dat b/tests/regression_tests/mgxs_library_hdf5/results_true.dat
@@ -97,8 +97,8 @@ domain=1 type=chi-prompt
 [1.00000000e+00 0.00000000e+00]
 [1.03333203e-01 0.00000000e+00]
 domain=1 type=inverse-velocity
-[5.72461488e-08 3.00999716e-06]
-[2.80644406e-09 1.80993425e-07]
+[5.72461488e-08 3.00999757e-06]
+[2.80644407e-09 1.80993426e-07]
 domain=1 type=prompt-nu-fission
 [5.64594959e-03 1.32859920e-01]
 [3.68364598e-04 7.79430310e-03]

diff --git a/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat b/tests/regression_tests/mgxs_library_no_nuclides/results_true.dat
@@ -143,7 +143,7 @@ chi-prompt
 inverse-velocity
    material  group in nuclide          mean     std. dev.
 1         1         1   total  5.956290e-08  2.255751e-09
-0         1         2   total  2.886630e-06  7.418455e-08
+0         1         2   total  2.886630e-06  7.418452e-08
 prompt-nu-fission
    material  group in nuclide      mean  std. dev.
 1         1         1   total  0.018067   0.000817
@@ -449,7 +449,7 @@ chi-prompt
 inverse-velocity
    material  group in nuclide          mean     std. dev.
 1         2         1   total  6.076563e-08  2.727519e-09
-0         2         2   total  2.996176e-06  1.110140e-07
+0         2         2   total  2.996176e-06  1.110139e-07
 prompt-nu-fission
    material  group in nuclide  mean  std. dev.
 1         2         1   total   0.0        0.0

diff --git a/tests/regression_tests/mgxs_library_nuclides/results_true.dat b/tests/regression_tests/mgxs_library_nuclides/results_true.dat
@@ -1 +1 @@
-c4a4cb4e00f09ef62222a0e66df817af87bf324a2ac0e57a82ff1337535a223c7077d660a60d0ce9ac7113c47d37b3296347f3ba212558ff09ef5fc586e1dd28
+4ae3b5a70ad72b1be261aee3ab19e0261d1c12f4d4ca50b712d1ba76041bd5387c69fa9fa326619ad588db206cd9aaf464b0025d71ea9b8b1137af0102bca87f
diff --git a/tests/regression_tests/mgxs_library_specific_nuclides/results_true.dat b/tests/regression_tests/mgxs_library_specific_nuclides/results_true.dat
@@ -1 +1 @@
-3e86542d1166b8a0bcc61742b88c9e931d63733477f3274b32b4da64d8f05413fa6330d5615f44f24735c2c2f77d3071f3dce38faba284c321ceab81c1064480
+08d5c199c51496f86fdd739bf7ee0e143a9a159da0f4d364ec970557e5c1fc92a202d906dcae91812e665fd2e88dd7db1e4913ef6b91f456f23b52093c83f483
diff --git a/tests/regression_tests/tallies/results_true.dat b/tests/regression_tests/tallies/results_true.dat
@@ -1 +1 @@
-d1decdbec6cb59df91ba5c42cb37a04f413a34fa7faf7ad1eecfd7d53a14af8cb65b58335c4ede4e88f6d9ab35a1746251983cc991d74eab3e678887c84183bd
+ddfbb0a6f5498eb8ff33bb10beb64e244b015861919eb4837bd82855e5fd87c3ff97dfa382d3d60afa81c48d9f857fba8e0b99263e7b35587f8adb75be1fc0ec