Random Ray Transport

[jtramm]

Introduction of Basic Random Ray Solver Mode in OpenMC

What is Random Ray?

Random ray is a stochastic transport method, closely related to the deterministic Method of Characteristics (MOC). Rather than each ray representing a single neutron as in Monte Carlo, it represents a characteristic line through the reactor upon which the transport equation can be written as an ordinary differential equation that can be solved analytically (although with discretization required in energy space, making it a multigroup method). The behavior of the governing transport equation can be approximated by solving along many characteristic tracks (rays) through the reactor. Unlike particles in Monte Carlo, rays in random ray or MOC are not affected by the material characteristics of the simulated problem -- rays are selected so as to explore the full simulation problem with a statistically equal distribution in space and angle.

An example of the transport process, showing two rays being processed, is given below:

random_ray_short.mov

Where can I read more about this method?

The below papers have some good background info for those that are interested. The third one may make for a good starting point as it introduces the method in the context of MC and compares/contrasts the two methods.

• First paper on random ray

• Larger random ray code results, performance and accuracy on 2D/3D C5G7

• In-depth comparison to multigroup Monte Carlo

Why is it being added to OpenMC?

There are a few good reasons:

1. The random ray solver complements the capabilities of MC nicely. One area that MC struggles with is maintaining accuracy in regions of low physical particle flux. Random ray, on the other hand, has approximately even variance throughout the entire global simulation domain, such that areas with low neutron flux are no less well known that areas of high neutron flux. Absent weight windows in MC, random ray can be several orders of magnitude faster than multigroup Monte Carlo in classes of problems where areas with low physical neutron flux need to be resolved. While MC uncertainty can be greatly improved with variance reduction techniques, they add some user complexity, and weight windows can often be expensive to generate via MC transport alone (e.g., MAGIC). While not yet implemented in this PR, I also plan to add capability so that the random ray solver can be used to generate fast and high quality weight windows via FW-CADIS that can then be used to accelerate convergence in MC. Early work in the field has shown significant speedup in weight window generation and weight window quality with random ray and FW-CADIS as compared to MAGIC.

2. In practical implementation terms, random ray is mechanically very similar to how Monte Carlo works, in terms of the process of ray tracing on CSG geometry and handling stochastic convergence, etc. In the original 1972 paper by Askew that introduces MOC (which random ray is a variant of), he stated:

One of the features of the method proposed [MoC] is that ... the tracking process needed to perform this operation is common to the proposed method ... and to Monte Carlo methods. Thus a single tracking routine capable of recognizing a geometric arrangement could be utilized to service all types of solution, choice being made depending which was more appropriate to the problem size and required accuracy.

This prediction holds up -- the additional requirements needed in OpenMC to handle random ray transport turned out to be fairly small, only in the range of 800 lines of code for the first implementation I did. Additionally, most of the solver was able to be added in a neatly "silo'd" fashion that doesn't complicate the MC portion of OpenMC. The current PR has increased in size to around 2400 lines of code, but much of that is due to the new example, regression test, and other boilerplate interface code etc. Thus, for relatively few additional lines of code, users will have capabilities to more efficiently handle certain types of simulation problems that are more computationally challenging to solve with MC, and will definitely have a faster transport method for generating weight windows.

3. It amortizes the code complexity in OpenMC for representing multigroup cross sections. There is a significant amount of interface code, documentation, and complexity in allowing OpenMC to generate and use multigroup XS data in its MGMC mode. Random ray allows the same multigroup data to be used, making full reuse of these existing capabilities.

Why not make this a standalone fork of OpenMC instead of adding it in?

This might have also been a reasonable approach. It could be argued that the new random ray solver creates bloat and might make it harder to maintain the repo. However, as the new solver is very small, and is written in a silo'd approach, there is very minimal added code complexity for anyone that is only concerned with developing MC features. As only a handful of very small changes are made to the main body of the repo (with the rest sequestered in the random_ray source folder), my guess is most developers won't even be aware that the solver was added unless they go looking for the changes.

Another key argument in favor of the random ray inclusion directly into OpenMC is to (eventually) enable relatively automatic generation of good weight windows, which are a key necessity in many fusion neutronics simulation problems. Requiring that weight windows be generated in an external code and maintaining the interface between the codes is likely much more work (both for users and developers) as compared to integrating weight window generation in the main OpenMC branch.

Why not add random ray to OpenMOC instead?

I had considered doing this, but I found it was a great deal easier to get this added to OpenMC as compared to OpenMOC. The reason being is that OpenMC is both inherently 3D and inherently stochastic, whereas the 3D extruded geometry ray tracing approach in OpenMOC would raise a lot of challenges in random ray for doing on-the-fly ray tracing. Additionally, it would take significant work in OpenMOC to allow it to handle inactive vs. active batches, accumulate unknowns with uncertainties, report standard deviations, etc, which we get for free in OpenMC.

How do I use the random ray mode?

Generally the inputs are all just about identical to what you'd see if doing a multigroup MC run. I've added a new pincell example to the repo to demonstrate its basic usage. The main changes you'll see are when configuring the settings for a run:

# Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.energy_mode = "multi-group"
settings.batches = 600
settings.inactive = 300
settings.particles = 50
settings.solver_type = 'random ray'
settings.random_ray_distance_inactive = 40.0
settings.random_ray_distance_active = 400.0

• particles setting now controls the number of rays
• solver_type is a new flag that can be used to specify random ray mode (defaults to Monte Carlo)
• random_ray_distance_inactive is the inactive (dead zone) distance that each ray travels before beginning integration (so as to build an estimate of the ray's starting angular flux on-the-fly).
• random_ray_distance_active is the total active distance that each ray travels while performing integration.

(If the dead zone and active length ideas are new to you, this paper has more info on these random ray concepts.)

This is about it in terms of the interface changes that are required to use the random ray mode. Beyond this, the only other major change you might want to do when converting a MGMC input deck into a random ray input deck is that you may need to subdivide cells into smaller regions so as to reduce the error associated with the flat source approximation used in this solver for random ray. In the future, we may add in higher order sources that will likely alleviate some of these requirements. We are also planning to add in the ability to automatically overlay a Cartesian mesh over the global geometry to more easily subdivide cells into flat source regions, but that is left for a future PR. The animation towards the top of the repo gives an idea of the level of subdivision typically used in light water reactor problems.

How are tallies handled?

Most tallies, filters, and scores that you would expect to work with a multigroup solver like random ray should work. E.g., you can define 3D mesh tallies with energy filters and flux, fission, and nu-fission scores, etc. There are some restrictions though. For starters, it is assumed that all filter mesh boundaries will conform to physical surface boundaries (or lattice boundaries) in the simulation geometry. It is acceptable for multiple cells (FSRs) to be contained within a filter mesh cell (e.g., pincell-level or assembly-level tallies should work), but it is currently left as undefined behavior if a single simulation cell is able to score to multiple filter mesh cells. In the future, we plan to add the capability to fully support mesh tallies, but for now this restriction needs to be respected.

How is plotting handled?

Visualization of geometry via openmc --plot is handled without any modifications, as the random ray solver uses the same geometry as in MC. However, any voxel plots that are defined will also be output at the end of the simulation generating .vtk files that can be directly read/plotting with Paraview. The purpose of this secondary (post simulation) plotting round is to output flat source region (FSR), power, and multigroup flux spectrum data on the voxel grid along with material information.

Why is plotting done like this?

In fixed source MC, by default the only thing we know after a simulation is the escape fraction. In a k-eigenvalue solve, by default all we know is the eigenvalue. Spatial flux information is left totally up to the user to record, and often fine-grained spatial meshes are considered costly/unnecessary, so it makes no sense in MC mode to try to attempt to plot any spatial flux or power info by default. Conversely, in random ray, the solver functions by estimating the multigroup source and flux spectrums in every (fairly small) cell each iteration. Thus, in random ray, in both fixed source and eigenvalue simulations, we always finish the simulation with a flux estimate for all areas. As such, it is much more common in MOC and other deterministic codes to plot in situ commonly as we always have spatial flux information available. Thus, OpenMC spits out .vtk files quickly at the end. In the future, all FSR data will be made available in the statepoint file, such that users will still have the ability to plot/manipulate it on the python end, although statepoint support is not yet available.

Why define a new plotter from scratch? Why not re-use the existing OpenMC plotter?

This is a fair point -- as the PR is being reviewed, I may try to see if there is an elegant way of making use of the existing plotter. The main issue that caused me to write my own is that the current voxel plotting is focused only on plotting geometry, and is only run at the beginning of the simulation, whereas for random ray we want to plot the spatial flux/power data after the simulation. In the interest of getting something done quickly, I simply reused a binary VTK plotter that I'd written previously, and found that it worked well and got the job done in ~150 lines of code.

What are the current limitations to the solver?

There are a bunch, currently! This is just the first basic implementation, more features and support will be added in subsequent PRs. However, current limitations are:

1. Mesh tallies currently have some limitations to them, as described earlier.

2. Only the following scores are currently supported: flux, total, fission, nu_fission, and events.

3. Only the following spatial tally filters are supported: cell, cell instance, distribcell, energy, material, mesh, and universe.

4. MGXS data must be isotropic and isothermal.

5. Fixed source transport is not yet supported, but will be soon.

6. You must always specify a single uniform SpatialBox isotropic source term for random ray, so that the solver knows where to sample ray starting points/angles from. The source term must not be biased, e.g., it should not be limited to fissionable regions, it should not be a point source. It should fill the entire spatial simulation domain.

7. While only voxel plots are accepted for generating outputs at the end of the simulation, all normal pre-simulation plots in plots.xml are supported.

8. MPI via domain replication is supported. However, domain replication does not scale well in random ray due to the need to all reduce the full scalar flux vector after each power iteration. Unlike in MC, where ranks only need to exchange a few fission sites with their neighbors to even the load, random ray requires a full all reduction of the scalar flux which is not scalable past a node or two. Domain replication was added for now though as it can be useful for very small problems (e.g., 2D pincell) that have very few cells causing memory contention issues. In these cases, domain replication can greatly speedup the runtime, so it does have its uses. Longer term, we plan on considering addition of domain decomposition (for both MC and random ray) so as to make the random ray solver more scalable.

9. User tallies will be written to statepoint files (and tallies.out, if enabled) as in a normal MC simulation. However, statepoints do not contain any other random ray state information as of yet (e.g., no cell-wise flux information or iteration source information is stored in the statepoint). As such, restart mode is not yet supported, and FSR flux/power information is not output anywhere yet. Results therefore are left to the contents of any tallies and any .vtk voxel plots.

What will be added in subsequent PRs? Is there a roadmap?

This PR is generally a "minimum viable product", with the basics of random ray added. A more robust feature set will be added in subsequent PR's over the next few months. The strategy of release a smaller basic PR first is to allow the community to try the solver out and give feedback before the layout of the solver starts to solidify as more advanced features are rolled out. Additionally, the smaller PR size makes it easier to review than if the branch was held until being "feature complete." Nonetheless, the current PR does have enough features to accomplish useful work. For instance, 2D C5G7 has been run and validated:

The main features slated for the next few months are:

• Shannon entropy reporting for evaluating inactive batch convergence.
• Cylindrical lattice definition for making FSR subdivision of pincells much easier and more performant.
• Improvements in mesh tallying capabilities so as to fully support all mesh types, and remove restriction that tally meshes conform to surfaces in the geometry.
• Performance optimization. Vectorization over energy groups at the inner loop may give significant speedup. Additionally, access of MGXS data (a frequent operation) is currently directly using the MGMC interface, which is not highly efficiently for the random ray use case. In the future, I may plan on either optimizing the MGMC interface or copying the data into a more ideal data structure for random ray's access pattern.
• Fixed source transport.
• Automatic mesh overlay to allow for easy subdivision of FSRs. This is likely to be a key feature for making fusion simulation problems that use DAGMC work with random ray.
• Statepoint storage of all flux/source data for postprocessing and for restarting.
• Adjoint solver
• In-depth theory and methods documentation.
• Weight window generation via FW-CADIS
• Full python interface

A few longer term stretch goals for random ray are also:

• Domain decomposition (for MC as well). Currently only MPI domain replication is supported.
• Streamlining of MGXS generation/usage workflow.

What is the runtime performance like?

The implementation in the current PR is lacking a few optimizations, but runtime info for random ray in OpenMC can be found in this paper. I'll plan on appending some formal C5G7 accuracy and timing results for this specific PR as a comment in the next few days.

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)

[paulromano]

Really nice work @jtramm! It's honestly amazing how few lines of codes it takes to add a fundamentally new solver in OpenMC. I already went through and made a few small fixes, updates, rewording, etc., so you may want to look through to make sure your OK with my updates. Otherwise, please see the relatively minor line comments below:

[gridley]

Looking forward to using this!

[yardasol]

Hi @jtramm, this is a really good first PR for this feature. I have given most of my comments in-line.

[jtramm]

Hi @yardasol -- if you're happy with the changes I've made, please switch your review status from "requested changes" to "approved". Thanks!

[yardasol]

Hi @jtramm,
Sorry for the delay, I've been at a conference. I think your changes to me feedback looks good. I've taken a look at the rest of the docs and think they are quite thorough. Great work!

[paulromano]

Good to go; thanks again for your work on this and especially for your patience @jtramm!