Add setting to optionally disable atomic relaxation

docs/source/io_formats/settings.rst

@@ -7,6 +7,19 @@ Settings Specification -- settings.xml
 All simulation parameters and miscellaneous options are specified in the
 settings.xml file.
 
+-------------------------------
+``<atomic_relaxation>`` Element
+-------------------------------
+
+The ``<atomic_relaxation>`` element determines whether the atomic relaxation
+cascade, the X-ray fluorescence photons and Auger electrons emitted when an
+inner-shell vacancy is filled, is simulated following photoelectric and
+incoherent (Compton) scattering interactions. Disabling this can speed up
+photon transport calculations where the detailed secondary particle cascade is
+not of interest.
+
+  *Default*: true
+
 ---------------------
 ``<batches>`` Element
 ---------------------

docs/source/usersguide/settings.rst

@@ -604,6 +604,13 @@ transport::
 
   settings.photon_transport = True
 
+Atomic relaxation (the cascade of fluorescence photons and Auger electrons
+emitted when an inner-shell vacancy is filled) is enabled by default whenever
+photon transport is on. It can be disabled using the
+:attr:`Settings.atomic_relaxation` attribute::
+
+  settings.atomic_relaxation = False
+
 The way in which OpenMC handles secondary charged particles can be specified
 with the :attr:`Settings.electron_treatment` attribute. By default, the
 :ref:`thick-target bremsstrahlung <ttb>` (TTB) approximation is used to generate

include/openmc/settings.h

@@ -77,6 +77,7 @@ extern "C" bool output_summary;      //!< write summary.h5?
 extern bool output_tallies;          //!< write tallies.out?
 extern bool particle_restart_run;    //!< particle restart run?
 extern "C" bool photon_transport;    //!< photon transport turned on?
+extern bool atomic_relaxation;       //!< atomic relaxation enabled?
 extern "C" bool reduce_tallies;      //!< reduce tallies at end of batch?
 extern bool res_scat_on;             //!< use resonance upscattering method?
 extern "C" bool restart_run;         //!< restart run?

openmc/settings.py

@@ -41,6 +41,10 @@ class Settings:
 
     Attributes
     ----------
+    atomic_relaxation : bool
+        Whether to simulate the atomic relaxation cascade (fluorescence photons
+        and Auger electrons) following photoelectric and incoherent scattering
+        interactions.
     batches : int
         Number of batches to simulate
     confidence_intervals : bool
@@ -402,6 +406,7 @@ def __init__(self, **kwargs):
         self._confidence_intervals = None
         self._electron_treatment = None
         self._photon_transport = None
+        self._atomic_relaxation = None
         self._plot_seed = None
         self._ptables = None
         self._uniform_source_sampling = None
@@ -663,6 +668,15 @@ def electron_treatment(self, electron_treatment: str):
                        electron_treatment, ['led', 'ttb'])
         self._electron_treatment = electron_treatment
 
+    @property
+    def atomic_relaxation(self) -> bool:
+        return self._atomic_relaxation
+
+    @atomic_relaxation.setter
+    def atomic_relaxation(self, atomic_relaxation: bool):
+        cv.check_type('atomic relaxation', atomic_relaxation, bool)
+        self._atomic_relaxation = atomic_relaxation
+
     @property
     def ptables(self) -> bool:
         return self._ptables
@@ -1631,6 +1645,11 @@ def _create_electron_treatment_subelement(self, root):
             element = ET.SubElement(root, "electron_treatment")
             element.text = str(self._electron_treatment)
 
+    def _create_atomic_relaxation_subelement(self, root):
+        if self._atomic_relaxation is not None:
+            element = ET.SubElement(root, "atomic_relaxation")
+            element.text = str(self._atomic_relaxation).lower()
+
     def _create_photon_transport_subelement(self, root):
         if self._photon_transport is not None:
             element = ET.SubElement(root, "photon_transport")
@@ -2129,6 +2148,11 @@ def _electron_treatment_from_xml_element(self, root):
         if text is not None:
             self.electron_treatment = text
 
+    def _atomic_relaxation_from_xml_element(self, root):
+        text = get_text(root, 'atomic_relaxation')
+        if text is not None:
+            self.atomic_relaxation = text in ('true', '1')
+
     def _energy_mode_from_xml_element(self, root):
         text = get_text(root, 'energy_mode')
         if text is not None:
@@ -2478,6 +2502,7 @@ def to_xml_element(self, mesh_memo=None):
         self._create_collision_track_subelement(element)
         self._create_confidence_intervals(element)
         self._create_electron_treatment_subelement(element)
+        self._create_atomic_relaxation_subelement(element)
         self._create_energy_mode_subelement(element)
         self._create_max_order_subelement(element)
         self._create_photon_transport_subelement(element)
@@ -2594,6 +2619,7 @@ def from_xml_element(cls, elem, meshes=None):
         settings._collision_track_from_xml_element(elem)
         settings._confidence_intervals_from_xml_element(elem)
         settings._electron_treatment_from_xml_element(elem)
+        settings._atomic_relaxation_from_xml_element(elem)
         settings._energy_mode_from_xml_element(elem)
         settings._max_order_from_xml_element(elem)
         settings._photon_transport_from_xml_element(elem)

src/photon.cpp

@@ -154,8 +154,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
 
     hid_t tgroup = open_group(rgroup, designator.c_str());
 
-    // Read binding energy energy and number of electrons if atomic relaxation
-    // data is present
+    // Read binding energy if atomic relaxation data is present
     if (attribute_exists(tgroup, "binding_energy")) {
       has_atomic_relaxation_ = true;
       read_attribute(tgroup, "binding_energy", shell.binding_energy);
@@ -174,7 +173,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
         i);
     cross_section = tensor::where(xs > 0, tensor::log(xs), 0);
 
-    if (object_exists(tgroup, "transitions")) {
+    if (settings::atomic_relaxation && object_exists(tgroup, "transitions")) {
       // Determine dimensions of transitions
       dset = open_dataset(tgroup, "transitions");
       auto dims = object_shape(dset);
@@ -206,9 +205,8 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   // Check the maximum size of the atomic relaxation stack
   auto max_size = this->calc_max_stack_size();
   if (max_size > MAX_STACK_SIZE && mpi::master) {
-    warning(fmt::format(
-      "The subshell vacancy stack in atomic relaxation can grow up to {}, but "
-      "the stack size limit is set to {}.",
+    warning(fmt::format("The subshell vacancy stack in atomic relaxation can "
+                        "grow up to {}, but the stack size limit is set to {}.",
       max_size, MAX_STACK_SIZE));
   }
 
@@ -231,7 +229,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
 
   // Map Compton subshell data to atomic relaxation data by finding the
   // subshell with the equivalent binding energy
-  if (has_atomic_relaxation_) {
+  if (settings::atomic_relaxation && has_atomic_relaxation_) {
     auto is_close = [](double a, double b) {
       return std::abs(a - b) / a < FP_REL_PRECISION;
     };

src/physics.cpp

@@ -355,7 +355,8 @@ void sample_photon_reaction(Particle& p)
     // Allow electrons to fill orbital and produce Auger electrons and
     // fluorescent photons. Since Compton subshell data does not match atomic
     // relaxation data, use the mapping between the data to find the subshell
-    if (i_shell >= 0 && element.subshell_map_[i_shell] >= 0) {
+    if (settings::atomic_relaxation && i_shell >= 0 &&
+        element.subshell_map_[i_shell] >= 0) {
       element.atomic_relaxation(element.subshell_map_[i_shell], p);
     }
 
@@ -427,7 +428,9 @@ void sample_photon_reaction(Particle& p)
 
         // Allow electrons to fill orbital and produce auger electrons
         // and fluorescent photons
-        element.atomic_relaxation(i_shell, p);
+        if (settings::atomic_relaxation) {
+          element.atomic_relaxation(i_shell, p);
+        }
         p.event() = TallyEvent::ABSORB;
         p.event_mt() = 533 + shell.index_subshell;
         p.wgt() = 0.0;

src/settings.cpp

@@ -62,6 +62,7 @@ bool output_summary {true};
 bool output_tallies {true};
 bool particle_restart_run {false};
 bool photon_transport {false};
+bool atomic_relaxation {true};
 bool reduce_tallies {true};
 bool res_scat_on {false};
 bool restart_run {false};
@@ -607,6 +608,11 @@ void read_settings_xml(pugi::xml_node root)
     }
   }
 
+  // Check for atomic relaxation
+  if (check_for_node(root, "atomic_relaxation")) {
+    atomic_relaxation = get_node_value_bool(root, "atomic_relaxation");
+  }
+
   // Number of bins for logarithmic grid
   if (check_for_node(root, "log_grid_bins")) {
     n_log_bins = std::stoi(get_node_value(root, "log_grid_bins"));

tests/regression_tests/atomic_relaxation/inputs_true.dat

@@ -0,0 +1,35 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <material id="1">
+      <density value="11.35" units="g/cm3"/>
+      <nuclide name="Pb208" ao="1.0"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" region="-1" universe="1"/>
+    <surface id="1" type="sphere" boundary="reflective" coeffs="0.0 0.0 0.0 1000000000.0"/>
+  </geometry>
+  <settings>
+    <run_mode>fixed source</run_mode>
+    <particles>10000</particles>
+    <batches>1</batches>
+    <source type="independent" strength="1.0" particle="photon">
+      <energy type="discrete">
+        <parameters>1000000.0 1.0</parameters>
+      </energy>
+    </source>
+    <electron_treatment>led</electron_treatment>
+    <atomic_relaxation>false</atomic_relaxation>
+    <photon_transport>true</photon_transport>
+  </settings>
+  <tallies>
+    <filter id="1" type="particle">
+      <bins>photon electron</bins>
+    </filter>
+    <tally id="1">
+      <filters>1</filters>
+      <scores>flux heating</scores>
+    </tally>
+  </tallies>
+</model>

tests/regression_tests/atomic_relaxation/results_true.dat

@@ -0,0 +1,9 @@
+tally 1:
+1.956204E+00
+3.826732E+00
+7.918768E+04
+6.270688E+09
+0.000000E+00
+0.000000E+00
+9.208123E+05
+8.478953E+11

tests/regression_tests/atomic_relaxation/test.py

@@ -0,0 +1,41 @@
+import openmc
+import pytest
+
+from tests.testing_harness import PyAPITestHarness
+
+
+@pytest.fixture
+def model():
+    mat = openmc.Material()
+    mat.add_nuclide('Pb208', 1.0)
+    mat.set_density('g/cm3', 11.35)
+
+    sphere = openmc.Sphere(r=1.0e9, boundary_type='reflective')
+    inside_sphere = openmc.Cell(fill=mat, region=-sphere)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([inside_sphere])
+
+    # Isotropic point source of 1 MeV photons at the origin
+    model.settings.source = openmc.IndependentSource(
+        particle='photon',
+        energy=openmc.stats.delta_function(1.0e6)
+    )
+
+    # Fixed-source photon transport with atomic relaxation disabled
+    model.settings.particles = 10000
+    model.settings.batches = 1
+    model.settings.photon_transport = True
+    model.settings.electron_treatment = 'led'
+    model.settings.atomic_relaxation = False
+    model.settings.run_mode = 'fixed source'
+
+    tally = openmc.Tally()
+    tally.filters = [openmc.ParticleFilter(['photon', 'electron'])]
+    tally.scores = ['flux', 'heating']
+    model.tallies = [tally]
+    return model
+
+
+def test_atomic_relaxation(model):
+    harness = PyAPITestHarness('statepoint.1.h5', model=model)
+    harness.main()

tests/unit_tests/test_settings.py

@@ -59,6 +59,7 @@ def test_export_to_xml(run_in_tmpdir):
     s.log_grid_bins = 2000
     s.photon_transport = False
     s.electron_treatment = 'led'
+    s.atomic_relaxation = False
     s.write_initial_source = True
     s.weight_window_checkpoints = {'surface': True, 'collision': False}
     source_region_mesh = openmc.RegularMesh()
@@ -147,6 +148,7 @@ def test_export_to_xml(run_in_tmpdir):
     assert s.log_grid_bins == 2000
     assert not s.photon_transport
     assert s.electron_treatment == 'led'
+    assert not s.atomic_relaxation
     assert s.write_initial_source
     assert len(s.volume_calculations) == 1
     vol = s.volume_calculations[0]