Random Ray AutoMagic Setup

[jtramm]

Description

This PR allows the user to almost entirely automate the process of converting an existing continuous energy model into a random ray model. Previously, the user had to:

1. Manually generate/provide multigroup cross section data
2. Convert their existing CE materials to MG materials by editing their specifications to remove isotopics and add the macroscopic XS data set
3. Populate all required random ray settings, which involves some knowledge of the random ray solver

This PR adds two new high level functions to the python interface.

model.convert_to_multigroup() automatically makes the conversion from CE -> MG. This function will invoke OpenMC internally to generate multigroup cross sections via continuous energy OpenMC simulations and store them to file. It will also replace the original continuous energy XS data sets with multigroup ones that point to the new library file. If a library file already exists, it will just use that rather than generate new XS data.

model.convert_to_random_ray() automatically populates the needed settings for random ray based on analysis of the model's geometry. This works well provided that the geometry's bounding box is reasonable.

The new interface thus works like:

  # Define regular openmc continuous energy model as normal
  model = openmc.pwr_pin_cell()

  # 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])]

Closes #3012.

How are MGXS Automatically Generated?

Prior to this PR, it was always up to the user to generate and provide their own MGXS library when running in multigroup Monte Carlo or with random ray. The new functionality in this PR provides an option for taking this responsibility away from the user, with a few options for different levels of fidelity vs. speed/generality are given. There will likely still be a lot of use cases where the user will still want to generate or provide their own MGXS library (e.g., when simulating a specific benchmark like C5G7 that includes specific multigroup data). However, the new functionality greatly improves the accessibility of the multigroup transport methods in OpenMC and should serve as a great starting point.

The new model.convert_to_multigroup() function takes a couple of optional arguments:

    def convert_to_multigroup(self, method="material-wise",
                              groups=openmc.mgxs.EnergyGroups(openmc.mgxs.GROUP_STRUCTURES['CASMO-2']),
                              nparticles=2000, overwrite_mgxs_library=False,
                              mgxs_fname: str = "mgxs.h5") -> None:

to allow the user to adjust the group structure, generation technique, and fidelity. This PR includes the following three methods for generating the MGXS:

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 distance from source	- 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

Notably, all of these methods are making big approximations, but they should provide enough fidelity for generating weight windows. In the future, perhaps we can add higher-fidelity MC generation methods (e.g., including transport correction, or source region specific XS data), or perhaps even faster deterministic generation methods.

How is the random ray solver configured?

The function convert_to_random_ray() will setup the random ray solver with some guesses based on the physical geometry of the simulation. It will use a minimum inactive length of 30.0 cm, though this will increase to the size of the longest chord length through the simulation domain (e.g., the distance between the lower_left and upper_right corners of the geometry bounding box). Active length is then 5x longer than the inactive length. Number of rays is a crude guess based on the general problem size -- the intention is that the user will likely want to adjust this after their first run when the miss rate is reported. It also uses the bounding box to form the sampling domain for the ray starting points. This last assumption typically works, but may not be valid if there are curved vacuum/reflective boundaries within the bounding box. In those cases, the user will need to manually adjust the bounding box to add a constraint that the rays be sampled within the physical domain (e.g., the top level universe containing the terminal boundary conditions and all geometry within them).

These are not meant to be the strictly "optimal" parameters, but are just half decent starting guesses that the user can then adjust as needed in response to observed results.

Documentation

I have now included a "quick start" guide at the top of the random ray section that mentions the new workflow, as I think most users will likely want to use these automation methods to get started. I have also updated the variance reduction section to use the simplified workflow, which greatly simplifies getting up and running with FW-CADIS.

Testing

I've added some new tests to cover the new workflow. It is just taking the standard CE pincell and applying the steps above, with variations for the three MGXS generation methods.

Validation

I have so far only done informal testing on the pin cell model. The default behavior tends to produce MGXS that result in a MGMC eigenvalue within about 100 pcm of the CEMC run, which is expected given the lack of transport correction in the moderator and the assumption of isotropy, etc. The random ray solver gives results very close to the MGMC results. As such, things appear to be working cleanly.

I also tested out the full core model in examples.py and it seems to run cleanly, though I haven't converged it yet.

Misc changes

• I also made a few tweaks to the python statepoint and summary classes so as to allow for manually closing their .h5 file handles. This was in response to complaints about hdf5 file locking due to running multiple calls to openmc.run() and subsequent statepoint loads needed by some of the MGXS workflows.

• I am reporting more data regarding settings in the stdout printout at the end of the simulation. This is helpful as the user may not be aware of what settings are used when they call model.convert_to_random_ray(), and the printout should make this more obvious if a wacky choice for a parameter is made by the function.

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 15) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[jeinstei]

I'm not sure if these "magic numbers" will become an issue, but is there a model.max_extents parameter that could be used for this instead? I can't imagine anyone wanting a larger simulation space, but it moves the magic into a parameter, at least.

[jtramm]

In this particular instance, the box is not placing any constraint on a user geometry -- instead this is just setting up some arbitrary reflective boundaries for a separate infinite medium calculation (i.e., just simulating an infinite extent of material) that is being run internally by the function. For the _generate_infinite_medium_mgxs() method, if there are 5 different materials in the user's real geometry, we will run 5 different simulations in an infinite medium, one for each material type. Theoretically we don't even need a reflective box at all, but OpenMC complains when there is no boundary condition given anywhere. Even if a particle does reach the box limit it just reflects back, so the 100k value doesn't affect the accuracy of the simulation. A box of length 1 cm could also be used without bias, but for optically thin materials like air this would increase the runtime due to all the added ray tracing.

[jeinstei]

Interesting -- I appreciate the rundown!