Add user setting for free gas threshold (#3593)

docs/source/io_formats/settings.rst

@@ -178,6 +178,16 @@ history-based parallelism.
 
   *Default*: false
 
+--------------------------------
+``<free_gas_threshold>`` Element
+--------------------------------
+
+The ``<free_gas_threshold>`` element specifies the energy multiplier, expressed
+in units of :math:`kT`, that determines when the free gas scattering approach is
+used for elastic scattering. Values must be positive.
+
+  *Default*: 400.0
+
 -----------------------------------
 ``<generations_per_batch>`` Element
 -----------------------------------

include/openmc/physics.h

@@ -10,13 +10,6 @@
 
 namespace openmc {
 
-//==============================================================================
-// Constants
-//==============================================================================
-
-// Monoatomic ideal-gas scattering treatment threshold
-constexpr double FREE_GAS_THRESHOLD {400.0};
-
 //==============================================================================
 // Non-member functions
 //==============================================================================

include/openmc/settings.h

@@ -147,6 +147,8 @@ extern std::unordered_set<int>
   source_write_surf_id; //!< Surface ids where sources will be written
 extern double source_rejection_fraction; //!< Minimum fraction of source sites
                                          //!< that must be accepted
+extern double free_gas_threshold;        //!< Threshold multiplier for free gas
+                                         //!< scattering treatment
 
 extern int
   max_history_splits; //!< maximum number of particle splits for weight windows

openmc/settings.py

@@ -84,6 +84,10 @@ class Settings:
         history-based parallelism.
 
         .. versionadded:: 0.12
+    free_gas_threshold : float
+        Energy multiplier (in units of :math:`kT`) below which the free gas
+        scattering treatment is applied for elastic scattering. If not
+        specified, a value of 400.0 is used.
     generations_per_batch : int
         Number of generations per batch
     ifp_n_generation : int
@@ -376,6 +380,7 @@ class Settings:
         self._seed = None
         self._stride = None
         self._survival_biasing = None
+        self._free_gas_threshold = None
 
         # Shannon entropy mesh
         self._entropy_mesh = None
@@ -1255,6 +1260,17 @@ class Settings:
         cv.check_less_than('source_rejection_fraction', source_rejection_fraction, 1)
         self._source_rejection_fraction = source_rejection_fraction
 
+    @property
+    def free_gas_threshold(self) -> float | None:
+        return self._free_gas_threshold
+
+    @free_gas_threshold.setter
+    def free_gas_threshold(self, free_gas_threshold: float | None):
+        if free_gas_threshold is not None:
+            cv.check_type('free gas threshold', free_gas_threshold, Real)
+            cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0)
+        self._free_gas_threshold = free_gas_threshold
+
     def _create_run_mode_subelement(self, root):
         elem = ET.SubElement(root, "run_mode")
         elem.text = self._run_mode.value
@@ -1733,6 +1749,11 @@ class Settings:
             element = ET.SubElement(root, "source_rejection_fraction")
             element.text = str(self._source_rejection_fraction)
 
+    def _create_free_gas_threshold_subelement(self, root):
+        if self._free_gas_threshold is not None:
+            element = ET.SubElement(root, "free_gas_threshold")
+            element.text = str(self._free_gas_threshold)
+
     def _eigenvalue_from_xml_element(self, root):
         elem = root.find('eigenvalue')
         if elem is not None:
@@ -2171,6 +2192,11 @@ class Settings:
         if text is not None:
             self.source_rejection_fraction = float(text)
 
+    def _free_gas_threshold_from_xml_element(self, root):
+        text = get_text(root, 'free_gas_threshold')
+        if text is not None:
+            self.free_gas_threshold = float(text)
+
     def to_xml_element(self, mesh_memo=None):
         """Create a 'settings' element to be written to an XML file.
 
@@ -2241,6 +2267,7 @@ class Settings:
         self._create_random_ray_subelement(element, mesh_memo)
         self._create_use_decay_photons_subelement(element)
         self._create_source_rejection_fraction_subelement(element)
+        self._create_free_gas_threshold_subelement(element)
 
         # Clean the indentation in the file to be user-readable
         clean_indentation(element)
@@ -2352,6 +2379,7 @@ class Settings:
         settings._random_ray_from_xml_element(elem, meshes)
         settings._use_decay_photons_from_xml_element(elem)
         settings._source_rejection_fraction_from_xml_element(elem)
+        settings._free_gas_threshold_from_xml_element(elem)
 
         return settings
 

src/finalize.cpp

@@ -85,6 +85,7 @@ int openmc_finalize()
   settings::time_cutoff = {INFTY, INFTY, INFTY, INFTY};
   settings::entropy_on = false;
   settings::event_based = false;
+  settings::free_gas_threshold = 400.0;
   settings::gen_per_batch = 1;
   settings::legendre_to_tabular = true;
   settings::legendre_to_tabular_points = -1;

src/physics.cpp

@@ -851,7 +851,7 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
 
     // otherwise, use free gas model
   } else {
-    if (E >= FREE_GAS_THRESHOLD * kT && nuc.awr_ > 1.0) {
+    if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) {
       return {};
     } else {
       sampling_method = ResScatMethod::cxs;

src/settings.cpp

@@ -126,6 +126,7 @@ SolverType solver_type {SolverType::MONTE_CARLO};
 std::unordered_set<int> sourcepoint_batch;
 std::unordered_set<int> statepoint_batch;
 double source_rejection_fraction {0.05};
+double free_gas_threshold {400.0};
 std::unordered_set<int> source_write_surf_id;
 int64_t ssw_max_particles;
 int64_t ssw_max_files;
@@ -651,6 +652,10 @@ void read_settings_xml(pugi::xml_node root)
       std::stod(get_node_value(root, "source_rejection_fraction"));
   }
 
+  if (check_for_node(root, "free_gas_threshold")) {
+    free_gas_threshold = std::stod(get_node_value(root, "free_gas_threshold"));
+  }
+
   // Survival biasing
   if (check_for_node(root, "survival_biasing")) {
     survival_biasing = get_node_value_bool(root, "survival_biasing");

tests/unit_tests/test_settings.py

@@ -80,6 +80,7 @@ def test_export_to_xml(run_in_tmpdir):
     s.max_particle_events = 100
     s.max_secondaries = 1_000_000
     s.source_rejection_fraction = 0.01
+    s.free_gas_threshold = 800.0
 
     # Make sure exporting XML works
     s.export_to_xml()
@@ -170,3 +171,4 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.random_ray['sample_method'] == 'halton'
     assert s.max_secondaries == 1_000_000
     assert s.source_rejection_fraction == 0.01
+    assert s.free_gas_threshold == 800.0