Random Ray AutoMagic Setup (#3351)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

.github/workflows/ci.yml

@@ -81,6 +81,11 @@ jobs:
       RDMAV_FORK_SAFE: 1
 
     steps:
+      - name: Setup cmake
+        uses: jwlawson/actions-setup-cmake@v2
+        with:
+          cmake-version: '3.31'
+
       - name: Checkout repository
         uses: actions/checkout@v4
         with:

docs/source/io_formats/settings.rst

@@ -487,6 +487,14 @@ found in the :ref:`random ray user guide <random_ray>`.
 
         :type:
           The type of the domain. Can be ``material``, ``cell``, or ``universe``.
+  
+  :diagonal_stabilization_rho:
+    The rho factor for use with diagonal stabilization. This technique is
+    applied when negative diagonal (in-group) elements are detected in
+    the scattering matrix of input MGXS data, which is a common feature
+    of transport corrected MGXS data.
+
+    *Default*: 1.0
 
 ----------------------------------
 ``<resonance_scattering>`` Element

docs/source/pythonapi/deplete.rst

@@ -78,20 +78,18 @@ A minimal example for performing depletion would be:
     >>> import openmc.deplete
     >>> geometry = openmc.Geometry.from_xml()
     >>> settings = openmc.Settings.from_xml()
-    >>> model = openmc.model.Model(geometry, settings)
+    >>> model = openmc.Model(geometry, settings)
 
     # Representation of a depletion chain
     >>> chain_file = "chain_casl.xml"
-    >>> operator = openmc.deplete.CoupledOperator(
-    ...     model, chain_file)
+    >>> operator = openmc.deplete.CoupledOperator(model, chain_file)
 
     # Set up 5 time steps of one day each
     >>> dt = [24 * 60 * 60] * 5
     >>> power = 1e6  # constant power of 1 MW
 
     # Deplete using mid-point predictor-corrector
-    >>> cecm = openmc.deplete.CECMIntegrator(
-    ...     operator, dt, power)
+    >>> cecm = openmc.deplete.CECMIntegrator(operator, dt, power)
     >>> cecm.integrate()
 
 Internal Classes and Functions

docs/source/usersguide/random_ray.rst

@@ -11,6 +11,76 @@ active batches <usersguide_batches>`. However, there are a couple of settings
 that are unique to the random ray solver and a few areas that the random ray
 run strategy differs, both of which will be described in this section.
 
+.. _quick_start:
+
+-----------
+Quick Start
+-----------
+
+While this page contains a comprehensive guide to the random ray solver and
+its various parameters, the process of converting an existing continuous energy
+Monte Carlo model to a random ray model can be largely automated via convenience
+functions in OpenMC's Python interface::
+
+  # Define continuous energy model as normal
+  model = openmc.Model()
+  ...
+
+  # Convert model to multigroup (will auto-generate MGXS library if needed)
+  model.convert_to_multigroup()
+
+  # Convert model to random ray and initialize random ray parameters
+  # to reasonable defaults based on the specifics of the geometry
+  model.convert_to_random_ray()
+
+  # (Optional) Overlay source region decomposition mesh to improve fidelity of the
+  # random ray solver. Adjust 'n' for fidelity vs runtime.
+  n = 100
+  mesh = openmc.RegularMesh()
+  mesh.dimension = (n, n, n)
+  mesh.lower_left = model.geometry.bounding_box.lower_left
+  mesh.upper_right = model.geometry.bounding_box.upper_right
+  model.settings.random_ray['source_region_meshes'] = [(mesh, [model.geometry.root_universe])]
+
+  # (Optional) Improve fidelity of the random ray solver by enabling linear sources
+  model.settings.random_ray['source_shape'] = 'linear'
+
+  # (Optional) Increase the number of rays/batch, to reduce uncertainty
+  model.settings.particles = 500
+
+The above strategy first converts the continuous energy model to a multigroup
+one using the :meth:`openmc.Model.convert_to_multigroup` method. By default,
+this will internally run a coarsely converged continuous energy Monte Carlo
+simulation to produce an estimated multigroup macroscopic cross section set for
+each material specified in the model, and store this data into a multigroup
+cross section library file (``mgxs.h5``) that can be used by the random ray
+solver.
+
+The :meth:`openmc.Model.convert_to_random_ray` method enables random ray mode
+and performs an analysis of the model geometry to determine reasonable values
+for all required parameters. If default behavior is not satisfactory, the user
+can manually adjust the settings in the :attr:`~openmc.Settings.random_ray`
+dictionary in the :class:`openmc.Settings` as described in the sections below.
+
+Finally a few optional steps are shown. The first (recommended) step overlays a
+mesh over the geometry to create smaller source regions so that source
+resolution improves and the random ray solver becomes more accurate. Varying the
+mesh resolution can be used to trade off between accuracy and runtime.
+High-fidelity fission reactor simulation may require source region sizes below 1
+cm, while larger fixed source problems with some tolerance for error may be able
+to use source regions of 10 or 100 cm.
+
+We also enable linear sources, which can improve the accuracy of the random ray
+solver and/or allow for a much coarser mesh resolution to be overlaid. Finally,
+the number of rays per batch is adjusted. The goal here is to ensure that the
+source region miss rate is below 1%, which is reported by OpenMC at the end of
+the simulation (or before via a warning if it is very high).
+
+.. warning::
+    If using a mesh filter for tallying or weight window generation, ensure that
+    the same mesh is used for source region decomposition via
+    ``model.settings.random_ray['source_region_meshes']``.
+
 ------------------------
 Enabling Random Ray Mode
 ------------------------
@@ -557,29 +627,118 @@ variety of problem types (or through a multidimensional parameter sweep of
 design variables) with only modest errors and at greatly reduced cost as
 compared to using only continuous energy Monte Carlo.
 
+~~~~~~~~~~~~
+The Easy Way
+~~~~~~~~~~~~
+
+The easiest way to generate a multigroup cross section library is to use the
+:meth:`openmc.Model.convert_to_multigroup` method. This method will
+automatically output a multigroup cross section library file (``mgxs.h5``) from
+a continuous energy Monte Carlo model and alter the material definitions in the
+model to use these multigroup cross sections. An example is given below::
+
+  # Assume we already have a working continuous energy model
+  model.convert_to_multigroup(
+      method="material_wise",
+      groups="CASMO-2",
+      nparticles=2000,
+      overwrite_mgxs_library=False,
+      mgxs_path="mgxs.h5",
+      correction=None
+  )
+
+The most important parameter to set is the ``method`` parameter, which can be
+either "stochastic_slab", "material_wise", or "infinite_medium". An overview
+of these methods is given below:
+
+.. list-table:: Comparison of Automatic MGXS Generation Methods
+   :header-rows: 1
+   :widths: 10 30 30 30
+
+   * - Method
+     - Description
+     - Pros
+     - Cons
+   * - ``material_wise`` (default)
+     - * Higher Fidelity
+       * Runs a CE simulation with the original geometry and source, tallying
+         cross sections with a material filter.
+     - * Typically the most accurate of the three methods
+       * Accurately captures (averaged over the full problem domain)
+         both spatial and resonance self shielding effects
+     - * Potentially slower as the full geometry must be run
+       * If a material is only present far from the source and doesn't get tallied
+         to in the CE simulation, the MGXS will be zero for that material.
+   * - ``stochastic_slab``
+     - * Medium Fidelity
+       * Runs a CE simulation with a greatly simplified geometry, where materials
+         are randomly assigned to layers in a 1D "stochastic slab sandwich" geometry
+     - * Still captures resonant self shielding and resonance effects between materials
+       * Fast due to the simplified geometry
+       * Able to produce cross section data for all materials, regardless of how
+         far they are from the source in the original geometry
+     - * Does not capture most spatial self shielding effects, e.g., no lattice physics.
+   * - ``infinite_medium``
+     - * Lower Fidelity
+       * Runs one CE simulation per material independently. Each simulation is just
+         an infinite medium slowing down problem, with an assumed external source term.
+     - * Simple
+     - * Poor accuracy (no spatial information, no lattice physics, no resonance effects
+         between materials)
+       * May hang if a material has a k-infinity greater than 1.0
+
+When selecting a non-default energy group structure, you can manually define
+group boundaries or specify the name of a known group structure (a list of which
+can be found at :data:`openmc.mgxs.GROUP_STRUCTURES`). The ``nparticles``
+parameter can be adjusted upward to improve the fidelity of the generated cross
+section library. The ``correction`` parameter can be set to ``"P0"`` to enable
+P0 transport correction. The ``overwrite_mgxs_library`` parameter can be set to
+``True`` to overwrite an existing MGXS library file, or ``False`` to skip
+generation and use an existing library file.
+
+.. note::
+    MGXS transport correction (via setting the ``correction`` parameter in the
+    :meth:`openmc.Model.convert_to_multigroup` method to ``"P0"``) may
+    result in negative in-group scattering cross sections, which can cause
+    numerical instability. To mitigate this, during a random ray solve OpenMC
+    will automatically apply
+    `diagonal stabilization <https://doi.org/10.1016/j.anucene.2018.10.036>`_
+    with a :math:`\rho` default value of 1.0, which can be adjusted with the
+    ``settings.random_ray['diagonal_stabilization_rho']`` parameter.
+
+Ultimately, the methods described above are all just approximations.
+Approximations in the generated MGXS data will fundamentally limit the potential
+accuracy of the random ray solver. However, the methods described above are all
+useful in that they can provide a good starting point for a random ray
+simulation, and if more fidelity is needed the user may wish to follow the
+instructions below or experiment with transport correction techniques to improve
+the fidelity of the generated MGXS data.
+
+~~~~~~~~~~~~
+The Hard Way
+~~~~~~~~~~~~
+
 We give here a quick summary of how to produce a multigroup cross section data
 file (``mgxs.h5``) from a starting point of a typical continuous energy Monte
-Carlo input file. Notably, continuous energy input files define materials as a
-mixture of nuclides with different densities, whereas multigroup materials are
-simply defined by which name they correspond to in a ``mgxs.h5`` library file.
+Carlo model. Notably, continuous energy models define materials as a mixture of
+nuclides with different densities, whereas multigroup materials are simply
+defined by which name they correspond to in a ``mgxs.h5`` library file.
 
 To generate the cross section data, we begin with a continuous energy Monte
-Carlo input deck and add in the required tallies that will be needed to generate
-our library. In this example, we will specify material-wise cross sections and a
-two group energy decomposition::
+Carlo model and add in the tallies that are needed to generate our library. In
+this example, we will specify material-wise cross sections and a two-group
+energy decomposition::
 
   # Define geometry
   ...
-  ...
   geometry = openmc.Geometry()
   ...
-  ...
 
   # Initialize MGXS library with a finished OpenMC geometry object
   mgxs_lib = openmc.mgxs.Library(geometry)
 
   # Pick energy group structure
-  groups = openmc.mgxs.EnergyGroups(openmc.mgxs.GROUP_STRUCTURES['CASMO-2'])
+  groups = openmc.mgxs.EnergyGroups('CASMO-2')
   mgxs_lib.energy_groups = groups
 
   # Disable transport correction
@@ -587,7 +746,7 @@ two group energy decomposition::
 
   # Specify needed cross sections for random ray
   mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
-                      'nu-scatter matrix', 'multiplicity matrix', 'chi']
+                         'nu-scatter matrix', 'multiplicity matrix', 'chi']
 
   # Specify a "cell" domain type for the cross section tally filters
   mgxs_lib.domain_type = "material"
@@ -614,13 +773,13 @@ two group energy decomposition::
   ...
 
 When selecting an energy decomposition, you can manually define group boundaries
-or pick out a group structure already known to OpenMC (a list of which can be
-found at :class:`openmc.mgxs.GROUP_STRUCTURES`). Once the above input deck has
-been run, the resulting statepoint file will contain the needed flux and
-reaction rate tally data so that a MGXS library file can be generated. Below is
-the postprocessing script needed to generate the ``mgxs.h5`` library file given
-a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
-``summary.h5``) that resulted from running our previous example::
+or specify the name of known group structure (a list of which can be found at
+:data:`openmc.mgxs.GROUP_STRUCTURES`). Once the above model has been run, the
+resulting statepoint file will contain the needed flux and reaction rate tally
+data so that a MGXS library file can be generated. Below is the postprocessing
+script needed to generate the ``mgxs.h5`` library file given a statepoint file
+(e.g., ``statepoint.100.h5``) file and summary file (e.g., ``summary.h5``) that
+resulted from running our previous example::
 
   import openmc
 
@@ -628,10 +787,7 @@ a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
   geom = summary.geometry
   mats = summary.materials
 
-  statepoint_filename = 'statepoint.100.h5'
-  sp = openmc.StatePoint(statepoint_filename)
-
-  groups = openmc.mgxs.EnergyGroups(openmc.mgxs.GROUP_STRUCTURES['CASMO-2'])
+  groups = openmc.mgxs.EnergyGroups('CASMO-2')
   mgxs_lib = openmc.mgxs.Library(geom)
   mgxs_lib.energy_groups = groups
   mgxs_lib.correction = None
@@ -653,10 +809,10 @@ a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
   # Construct all tallies needed for the multi-group cross section library
   mgxs_lib.build_library()
 
-  mgxs_lib.load_from_statepoint(sp)
+  with openmc.StatePoint('statepoint.100.h5') as sp:
+      mgxs_lib.load_from_statepoint(sp)
 
-  names = []
-  for mat in mgxs_lib.domains: names.append(mat.name)
+  names = [mat.name for mat in mgxs_lib.domains]
 
   # Create a MGXS File which can then be written to disk
   mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=names)
@@ -665,8 +821,8 @@ a statepoint file (e.g., ``statepoint.100.h5``) file and summary file (e.g.,
   mgxs_file.export_to_hdf5("mgxs.h5")
 
 Notably, the postprocessing script needs to match the same
-:class:`openmc.mgxs.Library` settings that were used to generate the tallies,
-but otherwise is able to discern the rest of the simulation details from the
+:class:`openmc.mgxs.Library` settings that were used to generate the tallies but
+is otherwise  able to discern the rest of the simulation details from the
 statepoint and summary files. Once the postprocessing script is successfully
 run, the ``mgxs.h5`` file can be loaded by subsequent runs of OpenMC.
 
@@ -701,11 +857,11 @@ multigroup library instead of defining their isotopic contents, as::
     water_data = openmc.Macroscopic('Hot borated water')
 
     # Instantiate some Materials and register the appropriate Macroscopic objects
-    fuel= openmc.Material(name='UO2 (2.4%)')
+    fuel = openmc.Material(name='UO2 (2.4%)')
     fuel.set_density('macro', 1.0)
     fuel.add_macroscopic(fuel_data)
 
-    water= openmc.Material(name='Hot borated water')
+    water = openmc.Material(name='Hot borated water')
     water.set_density('macro', 1.0)
     water.add_macroscopic(water_data)
 

docs/source/usersguide/variance_reduction.rst

@@ -88,59 +88,66 @@ random ray mode can be found in the :ref:`Random Ray User Guide <random_ray>`.
     ray solver. A high level overview of the current workflow for generation of
     weight windows with FW-CADIS using random ray is given below.
 
-1. Produce approximate multigroup cross section data (stored in a ``mgxs.h5``
-   library). There is more information on generating multigroup cross sections
-   via OpenMC in the :ref:`multigroup materials <create_mgxs>` user guide, and a
-   specific example of generating cross section data for use with random ray in
-   the :ref:`random ray MGXS guide <mgxs_gen>`.
+1. Begin by making a deepy copy of your continuous energy Python model and then
+   convert the copy to be multigroup and use the random ray transport solver.
+   The conversion process can largely be automated as described in more detail
+   in the :ref:`random ray quick start guide <quick_start>`, summarized below::
 
-2. Make a copy of your continuous energy Python input file. You'll edit the new
-   file to work in multigroup mode with random ray for producing weight windows.
+    # Define continuous energy model
+    ce_model = openmc.pwr_pin_cell() # example, replace with your model
 
-3. Adjust the material definitions in your new multigroup Python file to utilize
-   the multigroup cross sections instead of nuclide-wise continuous energy data.
-   There is a specific example of making this conversion in the :ref:`random ray
-   MGXS guide <mgxs_gen>`.
+    # Make a copy to convert to multigroup and random ray
+    model = copy.deepcopy(ce_model)
 
-4. Configure OpenMC to run in random ray mode (by adding several standard random
-   ray input flags and settings to the :attr:`openmc.Settings.random_ray`
-   dictionary). More information can be found in the  :ref:`Random Ray User
-   Guide <random_ray>`.
+    # Convert model to multigroup (will auto-generate MGXS library if needed)
+    model.convert_to_multigroup()
 
-5. Add in a :class:`~openmc.WeightWindowGenerator` in a similar manner as for
+    # Convert model to random ray and initialize random ray parameters
+    # to reasonable defaults based on the specifics of the geometry
+    model.convert_to_random_ray()
+
+    # (Optional) Overlay source region decomposition mesh to improve fidelity of the
+    # random ray solver. Adjust 'n' for fidelity vs runtime.
+    n = 10
+    mesh = openmc.RegularMesh()
+    mesh.dimension = (n, n, n)
+    mesh.lower_left = model.geometry.bounding_box.lower_left
+    mesh.upper_right = model.geometry.bounding_box.upper_right
+    model.settings.random_ray['source_region_meshes'] = [(mesh, [model.geometry.root_universe])]
+
+    # (Optional) Improve fidelity of the random ray solver by enabling linear sources
+    model.settings.random_ray['source_shape'] = 'linear'
+
+    # (Optional) Increase the number of rays/batch, to reduce uncertainty
+    model.settings.particles = 500
+
+   If you need to improve the fidelity of the MGXS library, there is more
+   information on generating multigroup cross sections via OpenMC in the
+   :ref:`random ray MGXS guide <mgxs_gen>`.
+
+2. Add in a :class:`~openmc.WeightWindowGenerator` in a similar manner as for
    MAGIC generation with Monte Carlo and set the :attr:`method` attribute set to
    ``"fw_cadis"``::
 
-       # Define weight window spatial mesh
-       ww_mesh = openmc.RegularMesh()
-       ww_mesh.dimension = (10, 10, 10)
-       ww_mesh.lower_left = (0.0, 0.0, 0.0)
-       ww_mesh.upper_right = (100.0, 100.0, 100.0)
-
-       # Create weight window object and adjust parameters
+       # Create weight window object and adjust parameters, using the same mesh
+       # we used for source region decomposition
        wwg = openmc.WeightWindowGenerator(
            method='fw_cadis',
-           mesh=ww_mesh,
+           mesh=mesh,
            max_realizations=settings.batches
        )
 
        # Add generator to openmc.settings object
        settings.weight_window_generators = wwg
 
-
 .. warning::
     If using FW-CADIS weight window generation, ensure that the selected weight
     window mesh does not subdivide any source regions in the problem. This can
-    be ensured by assigning the weight window tally mesh to the root universe so
-    as to create source region boundaries that conform to the mesh, as in the
-    example below.
+    be ensured by using the same mesh for both source region subdivision (i.e.,
+    assigning to ``model.settings.random_ray['source_region_meshes']``) and for
+    weight window generation.
 
-::
-    
-    root = model.geometry.root_universe    
-    settings.random_ray['source_region_meshes'] = [(ww_mesh, [root])]
-
-6. When running your multigroup random ray input deck, OpenMC will automatically
+3. When running your multigroup random ray input deck, OpenMC will automatically
    run a forward solve followed by an adjoint solve, with a
    ``weight_windows.h5`` file generated at the end. The ``weight_windows.h5``
    file will contain FW-CADIS generated weight windows. This file can be used in

include/openmc/random_ray/flat_source_domain.h

@@ -58,6 +58,7 @@ public:
   SourceRegionHandle get_subdivided_source_region_handle(
     int64_t sr, int mesh_bin, Position r, double dist, Direction u);
   void finalize_discovered_source_regions();
+  void apply_transport_stabilization();
   int64_t n_source_regions() const
   {
     return source_regions_.n_source_regions();
@@ -71,6 +72,9 @@ public:
   // Static Data members
   static bool volume_normalized_flux_tallies_;
   static bool adjoint_; // If the user wants outputs based on the adjoint flux
+  static double
+    diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization
+                                 // for transport corrected MGXS data
 
   // Static variables to store source region meshes and domains
   static std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
@@ -127,6 +131,12 @@ public:
   std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
     source_region_map_;
 
+  // If transport corrected MGXS data is being used, there may be negative
+  // in-group scattering cross sections that can result in instability in MOC
+  // and random ray if used naively. This flag enables a stabilization
+  // technique.
+  bool is_transport_stabilization_needed_ {false};
+
 protected:
   //----------------------------------------------------------------------------
   // Methods

openmc/model/model.py

@@ -1,9 +1,12 @@
 from __future__ import annotations
 from collections.abc import Iterable, Sequence
+import copy
 from functools import lru_cache
 from pathlib import Path
 import math
 from numbers import Integral, Real
+import random
+import re
 from tempfile import NamedTemporaryFile, TemporaryDirectory
 import warnings
 
@@ -1439,3 +1442,482 @@ class Model:
             self.materials = openmc.Materials(
                 self.geometry.get_all_materials().values()
             )
+
+    def _generate_infinite_medium_mgxs(self, groups, nparticles, mgxs_path, correction):
+        """Generate a MGXS library by running multiple OpenMC simulations, each
+        representing an infinite medium simulation of a single isolated
+        material. A discrete source is used to sample particles, with an equal
+        strength spread across each of the energy groups. This is a highly naive
+        method that ignores all spatial self shielding effects and all resonance
+        shielding effects between materials.
+
+        Parameters
+        ----------
+        groups : openmc.mgxs.EnergyGroups
+            Energy group structure for the MGXS.
+        mgxs_path : path-like
+            Filename for the MGXS HDF5 file.
+        """
+        warnings.warn("The infinite medium method of generating MGXS may hang "
+                      "if a material has a k-infinity > 1.0.")
+        mgxs_sets = []
+        for material in self.materials:
+            openmc.reset_auto_ids()
+            model = openmc.Model()
+
+            # Set materials on the model
+            model.materials = [material]
+
+            # Settings
+            model.settings.batches = 100
+            model.settings.particles = nparticles
+            model.settings.run_mode = 'fixed source'
+
+            # Make a discrete source that is uniform over the bins of the group structure
+            n_groups = groups.num_groups
+            midpoints = []
+            strengths = []
+            for i in range(n_groups):
+                bounds = groups.get_group_bounds(i+1)
+                midpoints.append((bounds[0] + bounds[1]) / 2.0)
+                strengths.append(1.0)
+
+            energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
+            model.settings.source = openmc.IndependentSource(
+                space=openmc.stats.Point(), energy=energy_distribution)
+            model.settings.output = {'summary': True, 'tallies': False}
+
+            # Geometry
+            box = openmc.model.RectangularPrism(
+                100000.0, 100000.0, boundary_type='reflective')
+            name = material.name
+            infinite_cell = openmc.Cell(name=name, fill=material, region=-box)
+            infinite_universe = openmc.Universe(name=name, cells=[infinite_cell])
+            model.geometry.root_universe = infinite_universe
+
+            # Add MGXS Tallies
+
+            # Initialize MGXS library with a finished OpenMC geometry object
+            mgxs_lib = openmc.mgxs.Library(model.geometry)
+
+            # Pick energy group structure
+            mgxs_lib.energy_groups = groups
+
+            # Disable transport correction
+            mgxs_lib.correction = correction
+
+            # Specify needed cross sections for random ray
+            if correction == 'P0':
+                mgxs_lib.mgxs_types = [
+                    'nu-transport', 'absorption', 'nu-fission', 'fission',
+                    'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
+                ]
+            elif correction is None:
+                mgxs_lib.mgxs_types = [
+                    'total', 'absorption', 'nu-fission', 'fission',
+                    'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
+                ]
+
+            # Specify a "cell" domain type for the cross section tally filters
+            mgxs_lib.domain_type = "material"
+
+            # Specify the cell domains over which to compute multi-group cross sections
+            mgxs_lib.domains = model.geometry.get_all_materials().values()
+
+            # Do not compute cross sections on a nuclide-by-nuclide basis
+            mgxs_lib.by_nuclide = False
+
+            # Check the library - if no errors are raised, then the library is satisfactory.
+            mgxs_lib.check_library_for_openmc_mgxs()
+
+            # Construct all tallies needed for the multi-group cross section library
+            mgxs_lib.build_library()
+
+            # Create a "tallies.xml" file for the MGXS Library
+            mgxs_lib.add_to_tallies_file(model.tallies, merge=True)
+
+            # Run
+            statepoint_filename = model.run()
+
+            # Load MGXS
+            with openmc.StatePoint(statepoint_filename) as sp:
+                mgxs_lib.load_from_statepoint(sp)
+
+            # Create a MGXS File which can then be written to disk
+            mgxs_set = mgxs_lib.get_xsdata(domain=material, xsdata_name=name)
+            mgxs_sets.append(mgxs_set)
+
+        # Write the file to disk
+        mgxs_file = openmc.MGXSLibrary(energy_groups=groups)
+        for mgxs_set in mgxs_sets:
+            mgxs_file.add_xsdata(mgxs_set)
+        mgxs_file.export_to_hdf5(mgxs_path)
+
+    @staticmethod
+    def _create_stochastic_slab_geometry(materials, cell_thickness=1.0, num_repeats=100):
+        """Create a geometry representing a stochastic "sandwich" of materials in a
+        layered slab geometry. To reduce the impact of the order of materials in
+        the slab, the materials are applied to 'num_repeats' different randomly
+        positioned layers of 'cell_thickness' each.
+
+        Parameters
+        ----------
+        materials : list of openmc.Material
+            List of materials to assign. Each material will appear exactly num_repeats times,
+            then the ordering is randomly shuffled.
+        cell_thickness : float, optional
+            Thickness of each lattice cell in x (default 1.0 cm).
+        num_repeats : int, optional
+            Number of repeats for each material (default 100).
+
+        Returns
+        -------
+        geometry : openmc.Geometry
+            The constructed geometry.
+        box : openmc.stats.Box
+            A spatial sampling distribution covering the full slab domain.
+        """
+        if not materials:
+            raise ValueError("At least one material must be provided.")
+
+        num_materials = len(materials)
+        total_cells = num_materials * num_repeats
+        total_width = total_cells * cell_thickness
+
+        # Generate an infinite cell/universe for each material
+        universes = []
+        for i in range(num_materials):
+            cell = openmc.Cell(fill=materials[i])
+            universes.append(openmc.Universe(cells=[cell]))
+
+        # Make a list of randomized material idx assignments for the stochastic slab
+        assignments = list(range(num_materials)) * num_repeats
+        random.seed(42)
+        random.shuffle(assignments)
+
+        # Create a list of the (randomized) universe assignments to be used
+        # when defining the problem lattice.
+        lattice_entries = [universes[m] for m in assignments]
+
+        # Create the RectLattice for the 1D material variation in x.
+        lattice = openmc.RectLattice()
+        lattice.pitch = (cell_thickness, total_width, total_width)
+        lattice.lower_left = (0.0, 0.0, 0.0)
+        lattice.universes = [[lattice_entries]]
+        lattice.outer = universes[0]
+
+        # Define the six outer surfaces with reflective boundary conditions
+        rpp = openmc.model.RectangularParallelepiped(
+            0.0, total_width, 0.0, total_width, 0.0, total_width,
+            boundary_type='reflective'
+        )
+
+        # Create an outer cell that fills with the lattice.
+        outer_cell = openmc.Cell(fill=lattice, region=-rpp)
+
+        # Build the geometry
+        geometry = openmc.Geometry([outer_cell])
+
+        # Define the spatial distribution that covers the full cubic domain
+        box = openmc.stats.Box(*outer_cell.bounding_box)
+
+        return geometry, box
+
+    def _generate_stochastic_slab_mgxs(self, groups, nparticles, mgxs_path, correction) -> None:
+        """Generate MGXS assuming a stochastic "sandwich" of materials in a layered
+        slab geometry. While geometry-specific spatial shielding effects are not
+        captured, this method can be useful when the geometry has materials only
+        found far from the source region that the "material_wise" method would
+        not be capable of generating cross sections for. Conversely, this method
+        will generate cross sections for all materials in the problem regardless
+        of type. If this is a fixed source problem, a discrete source is used to
+        sample particles, with an equal strength spread across each of the
+        energy groups.
+
+        Parameters
+        ----------
+        groups : openmc.mgxs.EnergyGroups
+            Energy group structure for the MGXS.
+        mgxs_path : path-like
+            Filename for the MGXS HDF5 file.
+        """
+        openmc.reset_auto_ids()
+        model = openmc.Model()
+        model.materials = self.materials
+
+        # Settings
+        model.settings.batches = 200
+        model.settings.inactive = 100
+        model.settings.particles = nparticles
+        model.settings.output = {'summary': True, 'tallies': False}
+        model.settings.run_mode = self.settings.run_mode
+
+        # Stochastic slab geometry
+        model.geometry, spatial_distribution = Model._create_stochastic_slab_geometry(
+            model.materials)
+
+        # Make a discrete source that is uniform over the bins of the group structure
+        n_groups = groups.num_groups
+        midpoints = []
+        strengths = []
+        for i in range(n_groups):
+            bounds = groups.get_group_bounds(i+1)
+            midpoints.append((bounds[0] + bounds[1]) / 2.0)
+            strengths.append(1.0)
+
+        energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
+        model.settings.source = [openmc.IndependentSource(
+            space=spatial_distribution, energy=energy_distribution, strength=1.0)]
+
+        model.settings.output = {'summary': True, 'tallies': False}
+
+        # Add MGXS Tallies
+
+        # Initialize MGXS library with a finished OpenMC geometry object
+        mgxs_lib = openmc.mgxs.Library(model.geometry)
+
+        # Pick energy group structure
+        mgxs_lib.energy_groups = groups
+
+        # Disable transport correction
+        mgxs_lib.correction = correction
+
+       # Specify needed cross sections for random ray
+        if correction == 'P0':
+            mgxs_lib.mgxs_types = ['nu-transport', 'absorption', 'nu-fission', 'fission',
+                                   'consistent nu-scatter matrix', 'multiplicity matrix', 'chi']
+        elif correction is None:
+            mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission',
+                                   'consistent nu-scatter matrix', 'multiplicity matrix', 'chi']
+
+        # Specify a "cell" domain type for the cross section tally filters
+        mgxs_lib.domain_type = "material"
+
+        # Specify the cell domains over which to compute multi-group cross sections
+        mgxs_lib.domains = model.geometry.get_all_materials().values()
+
+        # Do not compute cross sections on a nuclide-by-nuclide basis
+        mgxs_lib.by_nuclide = False
+
+        # Check the library - if no errors are raised, then the library is satisfactory.
+        mgxs_lib.check_library_for_openmc_mgxs()
+
+        # Construct all tallies needed for the multi-group cross section library
+        mgxs_lib.build_library()
+
+        # Create a "tallies.xml" file for the MGXS Library
+        mgxs_lib.add_to_tallies_file(model.tallies, merge=True)
+
+        # Run
+        statepoint_filename = model.run()
+
+        # Load MGXS
+        with openmc.StatePoint(statepoint_filename) as sp:
+            mgxs_lib.load_from_statepoint(sp)
+
+        names = [mat.name for mat in mgxs_lib.domains]
+
+        # Create a MGXS File which can then be written to disk
+        mgxs_file = mgxs_lib.create_mg_library(xs_type='macro', xsdata_names=names)
+        mgxs_file.export_to_hdf5(mgxs_path)
+
+    def _generate_material_wise_mgxs(self, groups, nparticles, mgxs_path, correction) -> None:
+        """Generate a material-wise MGXS library for the model by running the
+        original continuous energy OpenMC simulation of the full material
+        geometry and source, and tally MGXS data for each material. This method
+        accurately conserves reaction rates totaled over the entire simulation
+        domain. However, when the geometry has materials only found far from the
+        source region, it is possible the Monte Carlo solver may not be able to
+        score any tallies to these material types, thus resulting in zero cross
+        section values for these materials. For such cases, the "stochastic
+        slab" method may be more appropriate.
+
+        Parameters
+        ----------
+        groups : openmc.mgxs.EnergyGroups
+            Energy group structure for the MGXS.
+        mgxs_path : str
+            Filename for the MGXS HDF5 file.
+        """
+        openmc.reset_auto_ids()
+        model = copy.deepcopy(self)
+        model.tallies = openmc.Tallies()
+
+        # Settings
+        model.settings.batches = 200
+        model.settings.inactive = 100
+        model.settings.particles = nparticles
+        model.settings.output = {'summary': True, 'tallies': False}
+
+        # Add MGXS Tallies
+
+        # Initialize MGXS library with a finished OpenMC geometry object
+        mgxs_lib = openmc.mgxs.Library(model.geometry)
+
+        # Pick energy group structure
+        mgxs_lib.energy_groups = groups
+
+        # Disable transport correction
+        mgxs_lib.correction = correction
+
+        # Specify needed cross sections for random ray
+        if correction == 'P0':
+            mgxs_lib.mgxs_types = [
+                'nu-transport', 'absorption', 'nu-fission', 'fission',
+                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
+            ]
+        elif correction is None:
+            mgxs_lib.mgxs_types = [
+                'total', 'absorption', 'nu-fission', 'fission',
+                'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'
+            ]
+
+        # Specify a "cell" domain type for the cross section tally filters
+        mgxs_lib.domain_type = "material"
+
+        # Specify the cell domains over which to compute multi-group cross sections
+        mgxs_lib.domains = model.geometry.get_all_materials().values()
+
+        # Do not compute cross sections on a nuclide-by-nuclide basis
+        mgxs_lib.by_nuclide = False
+
+        # Check the library - if no errors are raised, then the library is satisfactory.
+        mgxs_lib.check_library_for_openmc_mgxs()
+
+        # Construct all tallies needed for the multi-group cross section library
+        mgxs_lib.build_library()
+
+        # Create a "tallies.xml" file for the MGXS Library
+        mgxs_lib.add_to_tallies_file(model.tallies, merge=True)
+
+        # Run
+        statepoint_filename = model.run()
+
+        # Load MGXS
+        with openmc.StatePoint(statepoint_filename) as sp:
+            mgxs_lib.load_from_statepoint(sp)
+
+        names = [mat.name for mat in mgxs_lib.domains]
+
+        # Create a MGXS File which can then be written to disk
+        mgxs_file = mgxs_lib.create_mg_library(
+            xs_type='macro', xsdata_names=names)
+        mgxs_file.export_to_hdf5(mgxs_path)
+
+    def convert_to_multigroup(self, method="material_wise", groups='CASMO-2',
+                              nparticles=2000, overwrite_mgxs_library=False,
+                              mgxs_path: PathLike = "mgxs.h5", correction=None):
+        """Convert all materials from continuous energy to multigroup.
+
+        If no MGXS data library file is found, generate one using one or more
+        continuous energy Monte Carlo simulations.
+
+        Parameters
+        ----------
+        method : {"material_wise", "stochastic_slab", "infinite_medium"}, optional
+            Method to generate the MGXS.
+        groups : openmc.mgxs.EnergyGroups or str, optional
+            Energy group structure for the MGXS or the name of the group
+            structure (based on keys from openmc.mgxs.GROUP_STRUCTURES).
+        mgxs_path : str, optional
+            Filename of the mgxs.h5 library file.
+        correction : str, optional
+            Transport correction to apply to the MGXS. Options are None and
+            "P0".
+        """
+        if isinstance(groups, str):
+            groups = openmc.mgxs.EnergyGroups(groups)
+
+        # Make sure all materials have a name, and that the name is a valid HDF5
+        # dataset name
+        for material in self.materials:
+            if material.name is None:
+                material.name = f"material {material.id}"
+            material.name = re.sub(r'[^a-zA-Z0-9]', '_', material.name)
+
+        # If needed, generate the needed MGXS data library file
+        if not Path(mgxs_path).is_file() or overwrite_mgxs_library:
+            if method == "infinite_medium":
+                self._generate_infinite_medium_mgxs(
+                    groups, nparticles, mgxs_path, correction)
+            elif method == "material_wise":
+                self._generate_material_wise_mgxs(
+                    groups, nparticles, mgxs_path, correction)
+            elif method == "stochastic_slab":
+                self._generate_stochastic_slab_mgxs(
+                    groups, nparticles, mgxs_path, correction)
+            else:
+                raise ValueError(
+                    f'MGXS generation method "{method}" not recognized')
+        else:
+            print(f'Existing MGXS library file "{mgxs_path}" will be used')
+
+        # Convert all continuous energy materials to multigroup
+        self.materials.cross_sections = mgxs_path
+        for material in self.materials:
+            material.set_density('macro', 1.0)
+            material._nuclides = []
+            material._sab = []
+            material.add_macroscopic(material.name)
+
+        self.settings.energy_mode = 'multi-group'
+
+    def convert_to_random_ray(self):
+        """Convert a multigroup model to use random ray.
+
+        This method determines values for the needed settings and adds them to
+        the settings.random_ray dictionary so as to enable random ray mode. The
+        settings that are populated are:
+
+        - 'ray_source' (openmc.IndependentSource): Where random ray starting
+          points are sampled from.
+        - 'distance_inactive' (float): The "dead zone" distance at the beginning
+          of the ray.
+        - 'distance_active' (float): The "active" distance of the ray
+        - 'particles' (int): Number of rays to simulate
+
+        The method will determine reasonable defaults for each of the above
+        variables based on analysis of the model's geometry. The function will
+        have no effect if the random ray dictionary is already defined in the
+        model settings.
+        """
+        # If the random ray dictionary is already set, don't overwrite it
+        if self.settings.random_ray:
+            warnings.warn("Random ray conversion skipped as "
+                          "settings.random_ray dictionary is already set.")
+            return
+
+        if self.settings.energy_mode != 'multi-group':
+            raise ValueError(
+                "Random ray conversion failed: energy mode must be "
+                "'multi-group'. Use convert_to_multigroup() first."
+            )
+
+        # Helper function for detecting infinity
+        def _replace_infinity(value):
+            if np.isinf(value):
+                return 1.0 if value > 0 else -1.0
+            return value
+
+        # Get a bounding box for sampling rays. We can utilize the geometry's bounding box
+        # though for 2D problems we need to detect the infinities and replace them with an
+        # arbitrary finite value.
+        bounding_box = self.geometry.bounding_box
+        lower_left = [_replace_infinity(v) for v in bounding_box.lower_left]
+        upper_right = [_replace_infinity(v) for v in bounding_box.upper_right]
+        uniform_dist_ray = openmc.stats.Box(lower_left, upper_right)
+        rr_source = openmc.IndependentSource(space=uniform_dist_ray)
+        self.settings.random_ray['ray_source'] = rr_source
+
+        # For the dead zone and active length, a reasonable guess is the larger of either:
+        # 1) The maximum chord length through the geometry (as defined by its bounding box)
+        # 2) 30 cm
+        # Then, set the active length to be 5x longer than the dead zone length, for the sake of efficiency.
+        chord_length = np.array(upper_right) - np.array(lower_left)
+        max_length = max(np.linalg.norm(chord_length), 30.0)
+
+        self.settings.random_ray['distance_inactive'] = max_length
+        self.settings.random_ray['distance_active'] = 5 * max_length
+
+        # Take a wild guess as to how many rays are needed
+        self.settings.particles = 2 * int(max_length)

openmc/settings.py

@@ -187,6 +187,17 @@ class Settings:
             or openmc.Universe. The mesh will be applied to the listed domains
             to subdivide source regions so as to improve accuracy and/or conform
             with tally meshes.
+        :diagonal_stabilization_rho:
+            The rho factor for use with diagonal stabilization. This technique is
+            applied when negative diagonal (in-group) elements are detected in
+            the scattering matrix of input MGXS data, which is a common feature
+            of transport corrected MGXS data. The default is 1.0, which ensures
+            no negative diagonal elements are present in the iteration matrix and
+            thus stabilizes the simulation. A value of 0.0 will disable diagonal
+            stabilization. Values between 0.0 and 1.0 will apply a degree of
+            stabilization, which may be desirable as stronger diagonal stabilization
+            also tends to dampen the convergence rate of the solver, thus requiring
+            more iterations to converge.
 
         .. versionadded:: 0.15.0
     resonance_scattering : dict
@@ -1177,6 +1188,10 @@ class Settings:
             elif key == 'sample_method':
                 cv.check_value('sample method', value,
                                ('prng', 'halton'))
+            elif key == 'diagonal_stabilization_rho':
+                cv.check_type('diagonal stabilization rho', value, Real)
+                cv.check_greater_than('diagonal stabilization rho',
+                                      value, 0.0, True)
             else:
                 raise ValueError(f'Unable to set random ray to "{key}" which is '
                                  'unsupported by OpenMC')
@@ -2018,7 +2033,7 @@ class Settings:
         if elem is not None:
             self.random_ray = {}
             for child in elem:
-                if child.tag in ('distance_inactive', 'distance_active'):
+                if child.tag in ('distance_inactive', 'distance_active', 'diagonal_stabilization_rho'):
                     self.random_ray[child.tag] = float(child.text)
                 elif child.tag == 'source':
                     source = SourceBase.from_xml_element(child)

src/random_ray/flat_source_domain.cpp

@@ -30,6 +30,7 @@ RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ {
   RandomRayVolumeEstimator::HYBRID};
 bool FlatSourceDomain::volume_normalized_flux_tallies_ {false};
 bool FlatSourceDomain::adjoint_ {false};
+double FlatSourceDomain::diagonal_stabilization_rho_ {1.0};
 std::unordered_map<int, vector<std::pair<Source::DomainType, int>>>
   FlatSourceDomain::mesh_domain_map_;
 
@@ -1106,6 +1107,11 @@ void FlatSourceDomain::flatten_xs()
           double sigma_s =
             m.get_xs(MgxsType::NU_SCATTER, g_in, &g_out, NULL, NULL, t, a);
           sigma_s_.push_back(sigma_s);
+          // For transport corrected XS data, diagonal elements may be negative.
+          // In this case, set a flag to enable transport stabilization for the
+          // simulation.
+          if (g_out == g_in && sigma_s < 0.0)
+            is_transport_stabilization_needed_ = true;
         }
       } else {
         sigma_t_.push_back(0);
@@ -1427,4 +1433,53 @@ void FlatSourceDomain::finalize_discovered_source_regions()
   discovered_source_regions_.clear();
 }
 
+// This is the "diagonal stabilization" technique developed by Gunow et al. in:
+//
+// Geoffrey Gunow, Benoit Forget, Kord Smith, Stabilization of multi-group
+// neutron transport with transport-corrected cross-sections, Annals of Nuclear
+// Energy, Volume 126, 2019, Pages 211-219, ISSN 0306-4549,
+// https://doi.org/10.1016/j.anucene.2018.10.036.
+void FlatSourceDomain::apply_transport_stabilization()
+{
+  // Don't do anything if all in-group scattering
+  // cross sections are positive
+  if (!is_transport_stabilization_needed_) {
+    return;
+  }
+
+  // Apply the stabilization factor to all source elements
+#pragma omp parallel for
+  for (int64_t sr = 0; sr < n_source_regions(); sr++) {
+    int material = source_regions_.material(sr);
+    if (material == MATERIAL_VOID) {
+      continue;
+    }
+    for (int g = 0; g < negroups_; g++) {
+      // Only apply stabilization if the diagonal (in-group) scattering XS is
+      // negative
+      double sigma_s =
+        sigma_s_[material * negroups_ * negroups_ + g * negroups_ + g];
+      if (sigma_s < 0.0) {
+        double sigma_t = sigma_t_[material * negroups_ + g];
+        double phi_new = source_regions_.scalar_flux_new(sr, g);
+        double phi_old = source_regions_.scalar_flux_old(sr, g);
+
+        // Equation 18 in the above Gunow et al. 2019 paper. For a default
+        // rho of 1.0, this ensures there are no negative diagonal elements
+        // in the iteration matrix. A lesser rho could be used (or exposed
+        // as a user input parameter) to reduce the negative impact on
+        // convergence rate though would need to be experimentally tested to see
+        // if it doesn't become unstable. rho = 1.0 is good as it gives the
+        // highest assurance of stability, and the impacts on convergence rate
+        // are pretty mild.
+        double D = diagonal_stabilization_rho_ * sigma_s / sigma_t;
+
+        // Equation 16 in the above Gunow et al. 2019 paper
+        source_regions_.scalar_flux_new(sr, g) =
+          (phi_new - D * phi_old) / (1.0 - D);
+      }
+    }
+  }
+}
+
 } // namespace openmc

src/random_ray/random_ray_simulation.cpp

@@ -478,6 +478,9 @@ void RandomRaySimulation::simulate()
       // Add source to scalar flux, compute number of FSR hits
       int64_t n_hits = domain_->add_source_to_scalar_flux();
 
+      // Apply transport stabilization factors
+      domain_->apply_transport_stabilization();
+
       if (settings::run_mode == RunMode::EIGENVALUE) {
         // Compute random ray k-eff
         k_eff_ = domain_->compute_k_eff(k_eff_);
@@ -577,17 +580,23 @@ void RandomRaySimulation::print_results_random_ray(
     header("Simulation Statistics", 4);
     fmt::print(
       " Total Iterations                  = {}\n", settings::n_batches);
-    fmt::print(" Flat Source Regions (FSRs)        = {}\n", n_source_regions);
     fmt::print(
-      " FSRs Containing External Sources  = {}\n", n_external_source_regions);
+      " Number of Rays per Iteration      = {}\n", settings::n_particles);
+    fmt::print(" Inactive Distance                 = {} cm\n",
+      RandomRay::distance_inactive_);
+    fmt::print(" Active Distance                   = {} cm\n",
+      RandomRay::distance_active_);
+    fmt::print(" Source Regions (SRs)              = {}\n", n_source_regions);
+    fmt::print(
+      " SRs Containing External Sources   = {}\n", n_external_source_regions);
     fmt::print(" Total Geometric Intersections     = {:.4e}\n",
       static_cast<double>(total_geometric_intersections));
     fmt::print("   Avg per Iteration               = {:.4e}\n",
       static_cast<double>(total_geometric_intersections) / settings::n_batches);
-    fmt::print("   Avg per Iteration per FSR       = {:.2f}\n",
+    fmt::print("   Avg per Iteration per SR        = {:.2f}\n",
       static_cast<double>(total_geometric_intersections) /
         static_cast<double>(settings::n_batches) / n_source_regions);
-    fmt::print(" Avg FSR Miss Rate per Iteration   = {:.4f}%\n", avg_miss_rate);
+    fmt::print(" Avg SR Miss Rate per Iteration    = {:.4f}%\n", avg_miss_rate);
     fmt::print(" Energy Groups                     = {}\n", negroups);
     fmt::print(
       " Total Integrations                = {:.4e}\n", total_integrations);
@@ -613,6 +622,33 @@ void RandomRaySimulation::print_results_random_ray(
     std::string adjoint_true = (FlatSourceDomain::adjoint_) ? "ON" : "OFF";
     fmt::print(" Adjoint Flux Mode                 = {}\n", adjoint_true);
 
+    std::string shape;
+    switch (RandomRay::source_shape_) {
+    case RandomRaySourceShape::FLAT:
+      shape = "Flat";
+      break;
+    case RandomRaySourceShape::LINEAR:
+      shape = "Linear";
+      break;
+    case RandomRaySourceShape::LINEAR_XY:
+      shape = "Linear XY";
+      break;
+    default:
+      fatal_error("Invalid random ray source shape");
+    }
+    fmt::print(" Source Shape                      = {}\n", shape);
+    std::string sample_method =
+      (RandomRay::sample_method_ == RandomRaySampleMethod::PRNG) ? "PRNG"
+                                                                 : "Halton";
+    fmt::print(" Sample Method                     = {}\n", sample_method);
+
+    if (domain_->is_transport_stabilization_needed_) {
+      fmt::print(" Transport XS Stabilization Used   = YES (rho = {:.3f})\n",
+        FlatSourceDomain::diagonal_stabilization_rho_);
+    } else {
+      fmt::print(" Transport XS Stabilization Used   = NO\n");
+    }
+
     header("Timing Statistics", 4);
     show_time("Total time for initialization", time_initialize.elapsed());
     show_time("Reading cross sections", time_read_xs.elapsed(), 1);

src/settings.cpp

@@ -345,6 +345,15 @@ void get_run_parameters(pugi::xml_node node_base)
         }
       }
     }
+    if (check_for_node(random_ray_node, "diagonal_stabilization_rho")) {
+      FlatSourceDomain::diagonal_stabilization_rho_ = std::stod(
+        get_node_value(random_ray_node, "diagonal_stabilization_rho"));
+      if (FlatSourceDomain::diagonal_stabilization_rho_ < 0.0 ||
+          FlatSourceDomain::diagonal_stabilization_rho_ > 1.0) {
+        fatal_error("Random ray diagonal stabilization rho factor must be "
+                    "between 0 and 1");
+      }
+    }
   }
 }
 

tests/regression_tests/mgxs_library_ce_to_mg/test.py

@@ -40,8 +40,8 @@ class MGXSTestHarness(PyAPITestHarness):
 
         # Build MG Inputs
         # Get data needed to execute Library calculations.
-        sp = openmc.StatePoint(self._sp_name)
-        self.mgxs_lib.load_from_statepoint(sp)
+        with openmc.StatePoint(self._sp_name) as sp:
+            self.mgxs_lib.load_from_statepoint(sp)
         self._model.mgxs_file, self._model.materials, \
             self._model.geometry = self.mgxs_lib.create_mg_mode()
 
@@ -55,11 +55,6 @@ class MGXSTestHarness(PyAPITestHarness):
         self._model.export_to_model_xml()
         self._model.mgxs_file.export_to_hdf5()
 
-        # Enforce closing statepoint and summary files so HDF5
-        # does not throw an error during the next OpenMC execution
-        sp._f.close()
-        sp._summary._f.close()
-
         # Re-run MG mode.
         if config['mpi']:
             mpi_args = [config['mpiexec'], '-n', config['mpi_np']]

tests/regression_tests/mgxs_library_ce_to_mg_nuclides/test.py

@@ -40,8 +40,8 @@ class MGXSTestHarness(PyAPITestHarness):
 
         # Build MG Inputs
         # Get data needed to execute Library calculations.
-        sp = openmc.StatePoint(self._sp_name)
-        self.mgxs_lib.load_from_statepoint(sp)
+        with openmc.StatePoint(self._sp_name) as sp:
+            self.mgxs_lib.load_from_statepoint(sp)
         self._model.mgxs_file, self._model.materials, \
             self._model.geometry = self.mgxs_lib.create_mg_mode()
 
@@ -55,11 +55,6 @@ class MGXSTestHarness(PyAPITestHarness):
         self._model.export_to_model_xml()
         self._model.mgxs_file.export_to_hdf5()
 
-        # Enforce closing statepoint and summary files so HDF5
-        # does not throw an error during the next OpenMC execution
-        sp._f.close()
-        sp._summary._f.close()
-
         # Re-run MG mode.
         if config['mpi']:
             mpi_args = [config['mpiexec'], '-n', config['mpi_np']]

tests/regression_tests/random_ray_auto_convert/infinite_medium/inputs_true.dat

@@ -0,0 +1,64 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material depletable="true" id="1" name="UO2__2_4__">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="UO2__2_4__"/>
+    </material>
+    <material id="2" name="Zircaloy">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Zircaloy"/>
+    </material>
+    <material id="3" name="Hot_borated_water">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Hot_borated_water"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="Fuel" region="-1" universe="0"/>
+    <cell id="2" material="2" name="Cladding" region="1 -2" universe="0"/>
+    <cell id="3" material="3" name="Water" region="2 3 -4 5 -6" universe="0"/>
+    <surface coeffs="0 0 0.39218" id="1" name="Fuel OR" type="z-cylinder"/>
+    <surface coeffs="0 0 0.4572" id="2" name="Clad OR" type="z-cylinder"/>
+    <surface boundary="reflective" coeffs="-0.63" id="3" name="left" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="4" name="right" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-0.63" id="5" name="bottom" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="6" name="top" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
+        </space>
+      </source>
+      <distance_inactive>30.0</distance_inactive>
+      <distance_active>150.0</distance_active>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="0" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+      <source_shape>linear</source_shape>
+    </random_ray>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <upper_right>0.63 0.63</upper_right>
+    </mesh>
+  </settings>
+  </model>

tests/regression_tests/random_ray_auto_convert/infinite_medium/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+7.796949E-01 1.055316E-02

tests/regression_tests/random_ray_auto_convert/material_wise/inputs_true.dat

@@ -0,0 +1,64 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material depletable="true" id="1" name="UO2__2_4__">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="UO2__2_4__"/>
+    </material>
+    <material id="2" name="Zircaloy">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Zircaloy"/>
+    </material>
+    <material id="3" name="Hot_borated_water">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Hot_borated_water"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="Fuel" region="-1" universe="0"/>
+    <cell id="2" material="2" name="Cladding" region="1 -2" universe="0"/>
+    <cell id="3" material="3" name="Water" region="2 3 -4 5 -6" universe="0"/>
+    <surface coeffs="0 0 0.39218" id="1" name="Fuel OR" type="z-cylinder"/>
+    <surface coeffs="0 0 0.4572" id="2" name="Clad OR" type="z-cylinder"/>
+    <surface boundary="reflective" coeffs="-0.63" id="3" name="left" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="4" name="right" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-0.63" id="5" name="bottom" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="6" name="top" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
+        </space>
+      </source>
+      <distance_inactive>30.0</distance_inactive>
+      <distance_active>150.0</distance_active>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="0" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+      <source_shape>linear</source_shape>
+    </random_ray>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <upper_right>0.63 0.63</upper_right>
+    </mesh>
+  </settings>
+  </model>

tests/regression_tests/random_ray_auto_convert/material_wise/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+7.375068E-01 7.015839E-03

tests/regression_tests/random_ray_auto_convert/stochastic_slab/inputs_true.dat

@@ -0,0 +1,64 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material depletable="true" id="1" name="UO2__2_4__">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="UO2__2_4__"/>
+    </material>
+    <material id="2" name="Zircaloy">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Zircaloy"/>
+    </material>
+    <material id="3" name="Hot_borated_water">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Hot_borated_water"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="Fuel" region="-1" universe="0"/>
+    <cell id="2" material="2" name="Cladding" region="1 -2" universe="0"/>
+    <cell id="3" material="3" name="Water" region="2 3 -4 5 -6" universe="0"/>
+    <surface coeffs="0 0 0.39218" id="1" name="Fuel OR" type="z-cylinder"/>
+    <surface coeffs="0 0 0.4572" id="2" name="Clad OR" type="z-cylinder"/>
+    <surface boundary="reflective" coeffs="-0.63" id="3" name="left" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="4" name="right" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-0.63" id="5" name="bottom" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="6" name="top" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>10</batches>
+    <inactive>5</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
+        </space>
+      </source>
+      <distance_inactive>30.0</distance_inactive>
+      <distance_active>150.0</distance_active>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="0" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+      <source_shape>linear</source_shape>
+    </random_ray>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <upper_right>0.63 0.63</upper_right>
+    </mesh>
+  </settings>
+  </model>

tests/regression_tests/random_ray_auto_convert/stochastic_slab/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+7.499679E-01 8.107614E-03

tests/regression_tests/random_ray_auto_convert/test.py

@@ -0,0 +1,54 @@
+import os
+
+import openmc
+from openmc.examples import pwr_pin_cell
+from openmc import RegularMesh
+from openmc.utility_funcs import change_directory
+import pytest
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+@pytest.mark.parametrize("method", ["material_wise", "stochastic_slab", "infinite_medium"])
+def test_random_ray_auto_convert(method):
+    with change_directory(method):
+        openmc.reset_auto_ids()
+
+        # Start with a normal continuous energy model
+        model = pwr_pin_cell()
+
+        # Convert to a multi-group model
+        model.convert_to_multigroup(
+            method=method, groups='CASMO-2', nparticles=30,
+            overwrite_mgxs_library=False, mgxs_path="mgxs.h5"
+        )
+
+        # Convert to a random ray model
+        model.convert_to_random_ray()
+
+        # Set the number of particles
+        model.settings.particles = 100
+
+        # Overlay a basic 2x2 mesh
+        n = 2
+        mesh = RegularMesh()
+        mesh.dimension = (n, n)
+        bbox = model.geometry.bounding_box
+        mesh.lower_left = (bbox.lower_left[0], bbox.lower_left[1])
+        mesh.upper_right = (bbox.upper_right[0], bbox.upper_right[1])
+        model.settings.random_ray['source_region_meshes'] = [
+            (mesh, [model.geometry.root_universe])]
+
+        # Set the source shape to linear
+        model.settings.random_ray['source_shape'] = 'linear'
+
+        harness = MGXSTestHarness('statepoint.10.h5', model)
+        harness.main()

tests/regression_tests/random_ray_diagonal_stabilization/inputs_true.dat

@@ -0,0 +1,65 @@
+<?xml version='1.0' encoding='utf-8'?>
+<model>
+  <materials>
+    <cross_sections>mgxs.h5</cross_sections>
+    <material depletable="true" id="1" name="UO2__2_4__">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="UO2__2_4__"/>
+    </material>
+    <material id="2" name="Zircaloy">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Zircaloy"/>
+    </material>
+    <material id="3" name="Hot_borated_water">
+      <density units="macro" value="1.0"/>
+      <macroscopic name="Hot_borated_water"/>
+    </material>
+  </materials>
+  <geometry>
+    <cell id="1" material="1" name="Fuel" region="-1" universe="0"/>
+    <cell id="2" material="2" name="Cladding" region="1 -2" universe="0"/>
+    <cell id="3" material="3" name="Water" region="2 3 -4 5 -6" universe="0"/>
+    <surface coeffs="0 0 0.39218" id="1" name="Fuel OR" type="z-cylinder"/>
+    <surface coeffs="0 0 0.4572" id="2" name="Clad OR" type="z-cylinder"/>
+    <surface boundary="reflective" coeffs="-0.63" id="3" name="left" type="x-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="4" name="right" type="x-plane"/>
+    <surface boundary="reflective" coeffs="-0.63" id="5" name="bottom" type="y-plane"/>
+    <surface boundary="reflective" coeffs="0.63" id="6" name="top" type="y-plane"/>
+  </geometry>
+  <settings>
+    <run_mode>eigenvalue</run_mode>
+    <particles>100</particles>
+    <batches>20</batches>
+    <inactive>15</inactive>
+    <source particle="neutron" strength="1.0" type="independent">
+      <space type="box">
+        <parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
+      </space>
+      <constraints>
+        <fissionable>true</fissionable>
+      </constraints>
+    </source>
+    <energy_mode>multi-group</energy_mode>
+    <random_ray>
+      <source particle="neutron" strength="1.0" type="independent">
+        <space type="box">
+          <parameters>-0.63 -0.63 -1.0 0.63 0.63 1.0</parameters>
+        </space>
+      </source>
+      <distance_inactive>30.0</distance_inactive>
+      <distance_active>150.0</distance_active>
+      <source_region_meshes>
+        <mesh id="1">
+          <domain id="0" type="universe"/>
+        </mesh>
+      </source_region_meshes>
+      <source_shape>linear</source_shape>
+      <diagonal_stabilization_rho>0.5</diagonal_stabilization_rho>
+    </random_ray>
+    <mesh id="1">
+      <dimension>2 2</dimension>
+      <lower_left>-0.63 -0.63</lower_left>
+      <upper_right>0.63 0.63</upper_right>
+    </mesh>
+  </settings>
+  </model>

tests/regression_tests/random_ray_diagonal_stabilization/results_true.dat

@@ -0,0 +1,2 @@
+k-combined:
+7.152917E-01 1.430362E-02

tests/regression_tests/random_ray_diagonal_stabilization/test.py

@@ -0,0 +1,61 @@
+import os
+
+from openmc.examples import pwr_pin_cell
+from openmc import RegularMesh
+
+from tests.testing_harness import TolerantPyAPITestHarness
+
+
+class MGXSTestHarness(TolerantPyAPITestHarness):
+    def _cleanup(self):
+        super()._cleanup()
+        f = 'mgxs.h5'
+        if os.path.exists(f):
+            os.remove(f)
+
+
+def test_random_ray_diagonal_stabilization():
+    # Start with a normal continuous energy model
+    model = pwr_pin_cell()
+
+    # Convert to a multi-group model, with 70 group XS
+    # and transport correction enabled. This will generate
+    # MGXS data with some negatives on the diagonal, in order
+    # to trigger diagonal correction.
+    model.convert_to_multigroup(
+        method='material_wise', groups='CASMO-70', nparticles=30,
+        overwrite_mgxs_library=True, mgxs_path="mgxs.h5", correction='P0'
+    )
+
+    # Convert to a random ray model
+    model.convert_to_random_ray()
+
+    # Set the number of particles
+    model.settings.particles = 100
+
+    # Overlay a basic 2x2 mesh
+    n = 2
+    mesh = RegularMesh()
+    mesh.dimension = (n, n)
+    bbox = model.geometry.bounding_box
+    mesh.lower_left = (bbox.lower_left[0], bbox.lower_left[1])
+    mesh.upper_right = (bbox.upper_right[0], bbox.upper_right[1])
+    model.settings.random_ray['source_region_meshes'] = [
+        (mesh, [model.geometry.root_universe])]
+
+    # Set the source shape to linear
+    model.settings.random_ray['source_shape'] = 'linear'
+
+    # Explicitly set the diagonal stabilization rho (default is otherwise 1.0).
+    # Note that if we set this to 0.0 (thus distabling stabilization), the
+    # problem should fail due to instability, so this is actually a good test
+    # problem.
+    model.settings.random_ray['diagonal_stabilization_rho'] = 0.5
+
+    # If rho was 0.0, the instability would cause failure after iteration 14,
+    # so we go a little past that.
+    model.settings.inactive = 15
+    model.settings.batches = 20
+
+    harness = MGXSTestHarness('statepoint.20.h5', model)
+    harness.main()

tools/ci/gha-install-dagmc.sh

@@ -3,7 +3,7 @@
 set -ex
 
 # MOAB Variables
-MOAB_BRANCH='Version5.1.0'
+MOAB_BRANCH='5.5.1'
 MOAB_REPO='https://bitbucket.org/fathomteam/moab/'
 MOAB_INSTALL_DIR=$HOME/MOAB/
 
@@ -19,7 +19,7 @@ cd $HOME
 mkdir MOAB && cd MOAB
 git clone -b $MOAB_BRANCH $MOAB_REPO
 mkdir build && cd build
-cmake ../moab -DENABLE_HDF5=ON -DENABLE_NETCDF=ON -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$MOAB_INSTALL_DIR -DENABLE_BLASLAPACK=OFF
+cmake ../moab -DENABLE_HDF5=ON -DENABLE_NETCDF=ON -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$MOAB_INSTALL_DIR
 make -j && make -j install
 rm -rf $HOME/MOAB/moab $HOME/MOAB/build