Development and validation of fully open-source R2S shutdown dose rate capabilities in OpenMC

We present the first fully open-source capabilities for shutdown dose rate (SDR) calculations of fusion energy facilities based on the Rigorous 2-Step (R2S) methodology. These capabilities have been implemented in the OpenMC Monte Carlo particle transport code, building on its existing capabilities while also leveraging new features that have been added to the code to support SDR calculations, such as decay photon source generation. Each of the individual physics components in the R2S workflow—neutron transport, activation, decay photon source generation, and photon transport—have been verified through code-to-code comparisons with MCNP6.2 and FISPACT-II 4.0. These comparisons generally demonstrate excellent agreement between codes for each of the physics components. The full cell-based R2S workflow was validated by performing a simulation of the first experimental campaign from the Frascati Neutron Generator (FNG) ITER dose rate benchmark problem from the Shielding INtegral Benchmark Archive and Database (SINBAD). For short cooling times, the dose calculated by OpenMC agrees with the experimental measurements within the stated experimental uncertainties. For longer cooling times, an overprediction of the shutdown dose was observed relative to experiment, which is consistent with previous studies in the literature. Altogether, these features constitute a combination of capabilities in a single, open-source codebase to provide the fusion community with a readily-accessible option for SDR calculations and a platform for rapidly analyzing the performance of fusion technology.


Introduction
Fusion energy systems present unique technological challenges that must be overcome to realize economical, gridready fusion energy.In particular, the design and prototyping of nuclear components for capturing the fusion power of highenergy neutrons must begin in earnest to meet the timelines proposed by both public and private entities [1][2][3][4].To successfully design and deploy blanket and first wall components, a powerful and accessible nuclear software ecosystem is necessary to expedite the learning process and shorten design iteration timescales from the current norm.The assessment of nuclear and radiological effects in a reactor design, known as neutronics, significantly impacts downstream requirements during design and exhibits a number of computational challenges that frequently result in serious analysis bottlenecks [5][6][7][8][9][10].
For example, high-energy neutrons generated from a burning deuterium-tritium (D-T) plasma are crucial for power extraction and tritium breeding but also irradiate and activate reactor materials.The resulting radionuclides from activated materials emit high-energy photons as they decay, producing a radiation dose after shutdown, known as the shutdown dose rate (SDR).An accurate assessment of the SDR is needed to ensure that routine operations and maintenance procedures meet radiation safety standards for licensing.SDR calculations represent a comprehensive use case for exercising fusion neutronics workflows as they are complex, computationally intensive, and span the space of fusion nuclear analysis endpoints, including the coupled physics of neutron transport, transmutation, and photon transport.
For these reasons, we consider an SDR calculation as the 'measuring stick' for the capabilities of a nuclear analysis toolchain in the design of fusion energy systems.There are two primary methods for performing SDR calculations: the Rigorous 2-Step (R2S) method [11] and the Direct 1-Step (D1S) method [12,13] (several recent articles have proposed new methods [14,15] that improve upon R2S and D1S but have not yet been widely adopted to our knowledge).The R2S and D1S methods each have their own advantages and limitations, but both have been used for large-scale SDR calculations in support of ITER, JET, and other fusion experiments [11,[16][17][18][19][20][21].A number of SDR workflows for these methods have been created by various institutions, a summary of which is shown in table 1.This table highlights the complexity of these existing workflows, which often involve numerous independent codes, some of which are patched or developed privately by other teams and not made publicly available.Typically these codes are linked together by input/output file generation scripts that are often not maintained by the developers of the codes being linked together.The user effort required to couple these diverse tools together into a workable process is error-prone and heavily constrained by export control, making the workflow incompatible with commercial production schedules.Furthermore, the D1S method requires custom code and data library patches that are not officially supported by MCNP [29] or other production tools.Ultimately, in every current implementation, the underlying codes and coupling frameworks are developed and maintained by disparate teams.This leads to duplicated efforts across many institutions, blurs lines between code maintenance responsibilities, and erects barriers to transparent, sustainable software development.
As motivation grows to embrace open science and opensource software for scientific research and discovery, there is increasing pressure within the world-wide fusion research community to follow suit and benefit from the many welldocumented advantages it brings [30,31].Specifically, the use of open-source tools for nuclear analysis and fusion technology development is gaining traction, particularly within academia and private industry.The open-source Monte Carlo particle transport code, OpenMC [32], is one such tool that has received significant attention from fusion research institutions and private companies.Over the last few years, OpenMC has added many of the key features necessary for analyzing fusion systems.However, there have been some features lacking that have slowed further adoption.This work presents the final building blocks for fusion nuclear analysis workflows and their implementation in OpenMC to support SDR calculations by the R2S method.Each physics module of the workflow is validated against state-of-the-art codes MCNP6.2 [33] or FISPACT-II 4.0 [34].Finally, a demonstration of this workflow as applied to the Frascatti Neutron Generator (FNG) ITER dose rate benchmark from the Shielding INtegral Benchmark Archive and Database (SINBAD) database [35] is presented.

R2S methodology and OpenMC implementation
The R2S method involves a neutron transport calculation to converge the neutron flux over a spatial discretization of the model (cell or mesh), activation of materials within these discrete regions based on this neutron flux, production of a photon source based on activated material compositions, and finally a photon transport calculation to estimate the SDR at the cooling times of interest [11,27].As illustrated in figure 1, the user inputs enter into the four physics modules represented by yellow rectangles, with intermediate data products in blue parallelograms.Existing R2S workflows employ a number of different codes for these physics modules as shown in table 1 and generate intermediate data products in a wide variety of formats.The novelty of our approach is that OpenMC serves as the solver for every physics module.The data products are a Other incarnations of R2S and D1S methods such as R2Smesh [27] and pyD1S [28] have been omitted as they are less widely used.
b GVR stands for global variance reduction.c Requires custom closed-source patch to source code.packaged in a standard format that can be processed easily via the OpenMC Python API.The physics solvers used in the R2S method-neutron transport, activation, photon source generation, and photon transport-are fully supported in OpenMC thanks to enhancements in recent releases [36,37].The additions presented herein further develop these capabilities and automate the R2S workflow so that a user can evaluate the SDR without manual information transfer from one solver to the next.The following subsections describe the implementation of the R2S workflow in OpenMC.

Neutron transport and activation
The first two steps of the R2S workflow, neutron transport and activation, are handled together by OpenMC's depletion module [38].The depletion module was originally written for simulating fuel depletion in a fission reactor.Fundamentally, the depletion module solves for the material compositions as a function of time resulting from the transmutation of nuclides by nuclear reactions and radioactive decay inside an irradiated environment, which can be expressed as where N i (t) ≡ density of nuclide i at time t σ i (E, t) ≡ transmutation cross section for nuclide i at energy E and time t ϕ (E, t) ≡ neutron flux at energy E and time t f j →i (E) ≡ fraction of transmutation reactions in nuclide j that produce nuclide i λ j →i ≡ decay constant for decay modes in nuclide j that produce nuclide i n ≡ total number of nuclides.
Equation (1) states that the rate of change of N i is equal to the production rate minus the loss rate.Because the equation for N i depends on the density of other nuclides, {N j | j ̸ = i}, equation (1) represents a system of first-order ordinary differential equations (ODEs).In order to form a proper initial value problem, the nuclide densities at time 0 must be specified: ( These equations can be written more compactly in matrix notation as where and A(n, t) ∈ R n×n is referred to as the burnup matrix containing the decay and transmutation coefficients.Note that the burnup matrix depends on n because the solution to the transport equation depends on the nuclide densities, thereby making equation (3) nonlinear.
The principal difficulty in solving equation (3) for a general case is that the elements of A are not constant with respect to time.Evaluating A for a given n and t is equivalent to solving the neutron transport equation (NTE) to obtain transmutation reaction rates and combining that information with decay data to construct the matrix of coefficients.Although the neutron transport equation is itself time-dependent, it is generally assumed that the timescale over which material compositions and the neutron source change is sufficiently long such that the transport equation can be solved as a steady state equation.For many fusion systems, the neutron source changes on timescales much longer than the neutron lifetime and thus this assumption is well satisfied.
When A is treated constant with respect to n and t, equation (3) admits a closed-form solution, n (t) = exp (At) n 0 . (5) Accordingly, solving equation (3) usually involves three separate components: (i) A method for solving the steady-state NTE to evaluate the matrix of coefficients, A(n, t), for a given n and t; (ii) A numerical method for integrating equation (3) forward in time, usually in a form that involves taking one or more matrix exponentials; and (iii) A method to evaluate the matrix exponential or, equivalently, the action of a matrix exponential on a vector.
In OpenMC, the method for solving the steady-state NTE is called the operator, and the method for numerically integrating forward in time is simply called the integrator.OpenMC has multiple integrators that can be selected by a user; full details can be found in [38].For all integrators, it is necessary to break up the total time of interest into discrete timesteps, each of which is assumed to correspond to a constant A(n, t).Time discretization is necessary for several reasons.First, the numerical methods themselves require discretization.Additionally, the source strength may be pulsed (piecewise constant), which places a constraint on the length of the timesteps since they cannot be longer than an individual pulse.Finally, the flux spectrum may change due to the change in material compositions over time, which means that a solution of the NTE at the beginning of irradiation may not be valid for the entire length of the irradiation.Each integrator also relies on the ability to evaluate the action of the matrix exponential on a vector.The method for this evaluation chosen by OpenMC is currently the incomplete partial fraction Chebyshev Rational Approximation Method (CRAM) [39].
As mentioned in [38], the time integration logic is completely separate from the logic for solving the NTE, so in principle the same integrators can be used with different operators.For the R2S workflow, we have taken advantage of this separation of concerns to provide multiple ways of evaluating the burnup matrix at each timestep; all methods still rely on OpenMC's underlying particle transport solver, but they differ in how the transmutation reaction rates are evaluated.The most accurate-and most expensive-approach is to solve the steady-state NTE at each timestep to determine fluxes and reaction rates.We refer to this as a fully coupled transportdepletion simulation, and in OpenMC it is enabled by the CoupledOperator Python class.For the R2S workflow, in most cases the change in material composition due to activation does not cause an appreciable change in the flux spectrum or reaction rates.Consequently, solving the NTE at every timestep is not necessary.Instead, a single neutron transport simulation can be carried out at the beginning of the workflow to determine fluxes and microscopic cross sections, which are then used to evaluate reaction rates at each timestep based on the source strength and the updated material compositions.So long as there is no appreciable 'burnup' of material that would induce changes in the flux spectra, this approach is much less costly while retaining accuracy.This approach in OpenMC is enabled by the IndependentOperator Python class.

Decay photon source generation
With predicted material compositions at each timestep, the next step in the R2S workflow is to compute the spatial and energy distribution of the photon source arising from the radioactive decay of any nuclides present in the activated material compositions.Let f i (E) be the number of photons emitted at energy E for the decay of a single atom of nuclide i.If λ i is the decay constant for nuclide i, the total number of photons emitted at energy E per unit time is then λ i f i (E).Accounting for all nuclides in the activated material composition, the total photon source is where V is the volume of the material.To calculate the decay photon source then, one needs access to decay constants and the decay photon distributions, f i (E).In the OpenMC implementation, these data are made available in the depletion chain file and are obtained from decay sublibrary files (MF = 8, MT = 457) within an ENDF-6 formatted library, such as ENDF/B or JEFF.Note that f i (E) may be represented as a series of discrete energies and intensities, a continuous distribution in E, or a combination thereof.
To calculate S γ (E) for a given activated material composition, OpenMC first combines any discrete and continuous photon spectra into a single distribution as well as combining the spectra for both gamma rays and x-rays since they are listed separately in a decay sublibrary evaluation, thus providing a single distribution for f i (E).Then, the individual distributions for each nuclide are combined into a single distribution weighted by the number of atoms of each nuclide as in equation (6).Finally, in the cell-based R2S workflow, a separate photon source is defined for each cell with an activated material composition.The source is defined to be spatially uniform over a bounding box covering the cell, and rejection is performed on sites within the volume so that the source is uniform only over the actual volume of the cell itself.
If M activated materials/cells are present in the problem, this results in M separate sources with intensities given by the corresponding integral of S γ (E).For a single source photon, it is necessary to sample among the M sources with probability proportional to their total intensity.OpenMC constructs a probability mass function over the M sources and then uses Walker alias sampling [40] to sample one of the M sources in O(1) time.Thus, even when M is very large, sampling the decay photon source does not have a negative impact on overall execution time.

Photon transport
With a decay photon source determined for every activated material at each cooling time, the last step in the R2S workflow is to carry out a photon transport simulation to determine the resulting dose in volumes of interest.All functionality needed for photon transport has already been implemented in OpenMC [41], so no additional development here was needed for the implementation of the R2S workflow.To determine the dose in a given volume, a normal flux tally is defined along with a cell filter and an energy-function filter that multiplies the tallied flux by an arbitrary function of energy at each tallying event (collision or tracklength).This allows any dose response function to be applied.

Verification of individual R2S physics modules
To assess the accuracy of the R2S workflow in OpenMC, two sets of comparisons have been carried out.First, to verify all the individual physics components and intermediate data presented in figure 1, we compare the outputs of OpenMC calculations against calculations from MCNP6.2 and FISPACT-II 4.0 (two of the most commonly used codes in other R2S code suites) given equivalent inputs.Second, in section 4, we compare OpenMC's R2S SDR evaluation on a benchmark problem against both a reference R2S code suite as well as experimental values.
For both verification of individual R2S physics modules as well as the full R2S workflow comparison in section 4, we have used the FNG ITER dose rate benchmark problem, hereafter referred to simply as the FNG dose benchmark.This benchmark is designed to emulate the irradiation of the ITER vacuum vessel.The experiment consists of several layers of water-equivalent PMMA (polymethyl methacrylate) embedded in a stainless steel block.An accelerated deuteron source strikes a tritium-titanium target, producing 14 MeV neutrons that irradiate the stainless steel block.After a pulsed irradiation was carried out over several days, the SDR was measured at various cooling times in a center cavity by both an active dosimeter and a thermoluminescent detector.This benchmark has been used extensively in past R2S workflow comparisons.A model of the geometry and materials for the FNG dose benchmark can be seen in figure 2. The majority of the model is comprised of steel and Perspex plates with a square cavity at the center of the assembly, which houses the activation foils.For detailed comparisons of neutron spectra and nuclide inventory, two cells are compared in future sections.These cells are highlighted in orange in figure 2. Cell 160 is close to the central cavity and is relatively unshielded from the neutron source.Cell 373 on the other hand is much further away from the central cavity and deep into the assembly.Further details on the specifics of the FNG dose benchmark can be found in the SINBAD database [35].

Neutron transport
The first step in the verification of the R2S workflow is neutron transport.There is already an extensive validation basis for the neutron transport capabilities in OpenMC, particularly for safety and criticality benchmarks [38,42,43].In addition, a number of recent publications have added to this body of validation data for fixed source neutron transport calculations [44,45].
To prepare an OpenMC model of the FNG dose benchmark, the existing MCNP model that is distributed as part of SINBAD was automatically converted using the MCNP model conversion utility for OpenMC [46].The conversion utility is able to seamlessly convert the geometry and material definitions but does not handle source specification or tallies.For the neutron transport verification, we have chosen to use an isotropic, monoenergetic 14.1 MeV point source at the origin.While this does not capture the true complexity of the angular and energy dependence of the neutron source, we consider it to be sufficient for the sake of an intercode comparison of resulting neutron flux spectra.
To obtain as close a comparison as possible, the material definitions in the MCNP model, which included several elemental compositions, were modified to exactly match the specification in OpenMC with natural elements expanded into individual isotopes.ENDF/B-VII.1 [47] neutron cross sections from the ENDF71x processed data library were used in both codes.The original ACE files from ENDF71x were directly converted to OpenMC's native HDF5 format, ensuring that there were no differences due to nuclear data processing (e.g.choice of resonance reconstruction tolerances in NJOY [48]).OpenMC was run with 10 7 particles per batch and 100 batches, whereas MCNP was run with 10 9 source particles to match the total number of particles simulated by OpenMC.In each code, the flux spectrum was tallied over the CCFE-709 bin energy group structure [34] in all 299 cells of the model.Furthermore, a set of weight windows generated by the FW-CADIS methodology [49] using the Attila deterministic solver within the Attila4MC [50] software were applied to both the MCNP and OpenMC neutron transport runs in order to achieve more uniform convergence of the flux spectra over the whole model and as a function of energy.
As shown in figure 3, we can see that provided the same model and underlying nuclear data libraries, the cell-averaged multigroup flux tally in OpenMC matches that computed by MCNP6.2 to an excellent degree.
Since it is difficult to convey the degree of agreement between the two transport solutions from the flux spectra of two individual cells, additional model-wide comparisons are provided for all 299 cells and all energy bins in figure 4. The top row of figure 4 (panels (a) and (b)) depicts the relative agreement between OpenMC and MCNP whereas the bottom row (panels (c) and (d)) depicts the absolute agreement.In figure 4(a) the probability distribution of the relative difference between OpenMC and MCNP for all cell-energy bin pairs with an MCNP relative error below 10% is provided.The reasoning behind this choice is simply that if the reference solution from MCNP is 'converged' by the metric of acceptable relative error in a particular energy bin we will compare that result to the result of OpenMC.Generally the energy bins for which the 10% relative error metric is exceeded correspond to thermal energies where there is very little flux.The same relative error filtering methodology is likewise applied to the data presented in figure 4(b).The cumulative distribution function (CDF) of the data in figure 4(a) is shown by the orange solid line in figure 4(b) and indicates that 68% of all cell-energy bins in the model converged by MCNP differ from the OpenMC tally by less than 1%, 97% differ by less than 10% and 99.5% differ by less than 20%.Whereas all the energy bins for all the cells are mixed into the histogram in figures 4(a) and (b) depicts similar information, but for individual cells.Each trajectory in figure 4(b) represents the CDF for the energy bins in a single cell as a function of the relative difference between OpenMC and MCNP.In this way we can see what fraction of the energy bins in a given cell exhibit a relative difference below a certain threshold.Including the CDF for each of the 299 cells in figure 4(b) allows us to state that >80% of the energy bins in all cells show less than a 10% difference between OpenMC and MCNP and >95.7% have differences less than 20%.
As the relative agreement between OpenMC and MCNP is given by figures 4(a) and (b), a sense of the convergence of both codes is given in figures 4(c) and (d). Figure 4(c) shows the relative error of the flux tally in each energy bin averaged over all cells for both MCNP (blue dashed line) and OpenMC (orange solid line), where the blue shaded region is bounded by the minimum and maximum relative error as a function of energy across all cells from the MCNP solution.Figure 4(d  shows the probability density function of the relative error in the flux tally for all cell-energy bin pairs in both OpenMC and MCNP.Given the excellent agreement between the MCNP and OpenMC flux tally results demonstrated by figures 4(a) and (b) and the distributions of the relative errors in figures 4(c) and (d), we can conclude that the transport solution computed by OpenMC agrees extremely well with the reference solution provided by MCNP and converges with similar characteristics as a function of energy and space over the entire model.With this confidence, we can proceed to use the flux spectra computed by OpenMC during this step as an input to FISPACT-II to perform our activation comparison.

Activation
Taking the cell-averaged flux spectra computed by OpenMC as exemplified by figure 3(a) for all non-void cells in the FNG dose benchmark we can begin our activation calculations with FISPACT-II 4.0.The irradiation scenario used in this verification is the campaign 1 schedule shown in table 2 and the nuclear data libraries used by FISPACT-II and the OpenMC depletion module were both ENDF/B-VII.1.For all non-void cells, a FISPACT-II calculation was run, including a cross-section collapse, decay data condense, library summary output, and the inventory calculation.To provide additional time resolution to the nuclide inventories during the activation calculations, each irradiation step and decay step of the campaign 1 schedule were broken down into 10 evenly spaced substeps purely for visualization purposes.
To obtain the most accurate comparison between FISPACT-II and OpenMC, a openmc.deplete.MicroXS object was generated directly from the FISPACT-II printlib.outfile for every cell.The printlib.outfile contains the one-group collapsed microscopic cross section values for the reaction pathways and nuclides present in the nuclear data library.The MicroXS object is simply the OpenMC Python representation of this mapping (i.e. a collection of cross sections, in barns, for a given set of nuclides and reactions).This data is then processed internally by OpenMC, along with an ENDF/B-VII.1 decay chain to produce the burnup matrix, A, from equation (3).This ensures both codes are using the same piecewise constant burnup matrix.By also using the same initial material composition, this implies that the only resulting differences in future compositions would be due to the different numerical integration techniques applied by each code.There are a number of differences between the two codes in this respect.OpenMC employs a wide range of numerical integration schemes (in this comparison, a simple predictor and by FISPACT-II in gray circular markers, showing excellent agreement for both long-and short-lived nuclides, as well as metastable states over a complex irradiation scenario. integrator was used) along with CRAM for evaluating the matrix exponential.FISPACT-II relies on the LSODES package, which directly integrates equation (3) with automated error control, avoiding the need to evaluate a matrix exponential.Nevertheless, the resulting inventory comparison is extremely good between the two methods as shown in figure 5.The jumps in 56 Mn activity around times of 1 d and 2 d are expected and correspond to the irradiation pulses listed in table 2.
In addition to the individual nuclide evolution in particular cells, the total activity of the FNG model was computed as a function of time.The comparison of activity levels between OpenMC and FISPACT-II can be seen in figure 6 where figure 6(a) depicts the absolute comparison and figure 6(b) depicts the relative difference between them.An increase in the total activity is observed during irradiation as dictated by the schedule in table 2. For the majority of time points and, importantly, all time points in the shutdown phase, the two agree to within a fraction of a percent.There are a few transients during the pulsed irradiation scenario where the difference reaches slightly above 1%.Some important details to note aside from the numerical integration where differences could arise between OpenMC and FISPACT-II include the handling of within-group reaction rates as well as complex breakup reactions.One substantial benefit available to transport coupled depletion codes like OpenMC is the ability to correctly compute (and preserve) within-group reaction rates.Computing the one-group microscopic cross sections and fluxes for populating the burnup matrix within a continuous-energy Monte Carlo (MC) transport code ensures that energy self-shielding of all nuclides within the same material is properly accounted for and that reaction rates are preserved.When the transport solver and depletion solver are only coupled through a multigroup flux spectrum, approximations must be made to correct for these effects.

Decay photon source generation
To verify the accuracy of the generation of the decay photon source used in the photon transport step, we have carried out several comparisons of decay photon spectra generated by OpenMC compared to spectra generated by FISPACT-II 4.0.First, we compared the decay photon spectra between OpenMC and FISPACT-II for a material representing 1 gram of a single nuclide and repeated this comparison for every unstable nuclide in the ENDF/B-VIII.0[51] decay sublibrary that had spectral data for photons (gamma rays and/or x-rays).OpenMC represents the decay photon source using the discrete lines, whereas FISPACT-II bins the discrete lines into a 24 energy group structure.Thus, to compare results, the spectra generated by OpenMC were binned onto the same energy group structure, after which we compared the relative difference in the energy production in each group in units of eV s −1 .Histogram of relative differences between OpenMC and FISPACT-II for the calculated average energy produced in each of the 24 energy groups from each unstable nuclide in ENDF/B-VIII.0 with a discrete spectrum.
Figure 7 shows a histogram of the relative differences over all energy groups for all nuclides.Note that in most cases, the binned decay photon spectra do not have non-zero values in all energy groups, and thus the total number of data points that were compared is less than 24 times the number of nuclides.The results in figure 7 demonstrate that the average energy produced in each group agrees very closely between OpenMC and FISPACT-II; in all cases the relative difference is less than 10 −5 .While these results do illustrate excellent agreement, it is worth noting possible sources of discrepancy.One known difference is that FISPACT-II estimates the contribution of high-energy decay electrons to the photon spectrum from bremsstrahlung radiation, whereas OpenMC currently does not account for bremsstrahlung production during the generation of a decay photon source.
Note that the comparison in figure 7 only includes nuclides for which the decay photon spectra were represented as discrete lines.In the case where the spectrum in an ENDF decay sublibrary evaluation is represented as a continuous distribution in energy, OpenMC will use the distribution as-is; however, FISPACT-II is unable to use the continuous distribution directly and instead relies on an approximation [34].For these cases, much larger relative differences would be expected.Out of the entire ENDF/B-VIII.0decay sublibrary, this occurred for 276 nuclides.Note that all of these nuclides are short-lived: 76% of them have half-lives of less than 1 s and the longest half-life observed over all 276 nuclides was 5.3 min.Thus, the impact of these differences on an SDR calculation is expected to be negligible.
In addition to the comparisons performed on the decay photon spectra from single nuclides, we also compared the resulting decay photon spectra from realistic, activated material compositions consisting of many nuclides to confirm that the photon spectra produced by OpenMC agrees with the FISPACT-II spectra when binned onto the same energy group structure.Figure 8 shows the discrete photon spectrum calculated by OpenMC for an activated steel composition from the FNG dose benchmark along with the same distribution binned onto the 24 energy group structure and compared to the spectrum generated by FISPACT-II.Again, excellent agreement is observed between OpenMC and FISPACT-II.

Photon transport and dose conversion
Finally, to verify the photon transport capabilities in OpenMC, we compare the dose calculated during a photon transport simulation of the FNG dose benchmark between OpenMC and MCNP6.2.For this comparison, the photon source was defined as a spatially uniform distribution over a single cell in the FNG dose benchmark with a binned energy distribution going up to 4 MeV.The dose in the central cavity was calculated in each code by tallying the flux multiplied by the ICRP-74 flux-to-dose conversion factors and binned on the CCFE-709 bin energy group structure.Each code simulated 10 8 .Both OpenMC and MCNP used ENDF/B-VIII.0photoatomic and atomic relaxation data; for MCNP, the eprdata14 library was used.
Figure 9 shows the dose in the central cavity as calculated by OpenMC and MCNP.One can observe that across the entire range of energies, the dose calculated by OpenMC matches the dose calculated by MCNP within uncertainty (two standard deviations).This comparison thus demonstrates that given the same photoatomic data, decay photon source, and flux-to-dose conversion factors, OpenMC can produce a resulting dose that is statistically identical with MCNP.

Validation of the OpenMC R2S workflow
This section presents a validation of the full cell-based R2S workflow in OpenMC against experimental measurements made for the FNG dose benchmark.Comparisons are also made to previous results using the MCR2S code system.The workflow was executed with roughly 200 lines of Python using no dependencies other than OpenMC and widely used third-party Python packages.We expect the number of lines required for a typical R2S application to decrease over time as the API is refined.Two different experimental campaigns were carried out at the FNG to support this benchmark.Simulation of the first experimental campaign, which was conducted between May and July 2000, was carried out with OpenMC's R2S capabilities to compare against the experimental measurements.
As with the comparisons in section 3, an OpenMC model of the FNG dose benchmark was produced through an automated MCNP model conversion utility for OpenMC.Note that separate MCNP models are provided in SINBAD for the irradiation period (for neutron transport) and the cooling period (for photon transport); accordingly two OpenMC models are used during the R2S calculation.The original MCNP model relies on a compiled source file to represent the neutron source distribution from D-T reactions.While it is possible to 'wrap' this source file and use it within OpenMC, we have chosen instead to generate a tabulated representation directly from the compiled source file that can be sampled in OpenMC using its standard distribution types for the sake of simplicity and reproducibility.The tabulated representation is openly available in a public repository containing associated research data [52].Note that the source is tabulated in energy and angle but the spatial distribution is left as a point source.The original compiled source includes a beam width of 0.7 cm in the x-z plane that is expected to have a very small effect on the overall results given that the full experiment extends about 1 m in each direction.
The first experimental campaign at FNG was carried out over two months.First, the stainless steel block was irradiated with the D-T neutron source for 18 h over three days.The full scenario is detailed in table 2. In the experiment, the dose rate after irradiation was measured continuously using a Geiger-Müller counter in the center cavity for 60 d.Additionally, reaction rates on 58 Ni were measured during irradiation using six Ni activation foils placed along the center cavity.
A cell-based R2S calculation of the campaign 1 experiment was carried out using OpenMC.All steps of the calculation were executed on 16 nodes of the Bebop cluster within Argonne's Laboratory Computing Resource Center.Each node on Bebop has two 18-core Intel Xeon E5-2695v4 processors.First, a neutron transport simulation with 10 8 source particles (100 batches of 10 6 particles) was carried out to determine one-group microscopic cross sections and neutron fluxes in each irradiated cell.Unique copies of the material composition were created for each of the irradiated cells and are then passed to the IndependentOperator class along with the irradiation schedule to carry out the activation calculation.To assess the dose rate, the activated composition was determined at the following cooling times: The R2S calculations were carried out with ENDF/B-VIII.0cross sections.A depletion chain was generated using a combination of the incident neutron, decay, and fission product yield ENDF sublibraries.The full depletion chain was then reduced to a special-purpose activation chain by starting from the nuclides that are present across all materials in the FNG dose benchmark and following their transmutation paths four times.This results in a depletion chain with 492 nuclides.A visual representation of the nuclides present in this reduced depletion chain is shown in figure 10 compared to a depletion chain with all nuclides from ENDF/B-VIII.0.
During the neutron transport step, tallies were included in the model for the (n,p) and (n,γ) reaction rates on 58 Ni for each of six activation foils.The calculated reaction rates on the Ni foils are presented in tables 3 and 4 along with the experimental values.To assess the impact of different nuclear data libraries, separate calculations have been performed with the ENDF/B-VIII.0and FENDL-3.2bdata libraries, each with a   total of 10 9 source neutrons.In the original MCNP neutron irradiation model, the cells intended for the six Ni foils contain air.To calculate the reaction rate, a tally was defined for the neutron flux multiplied by the 58 Ni (n,p) and (n,2n) microscopic cross sections.Thus, the resulting tally provides the 58 Ni reaction rates at a hypothetical concentration of 1 atom/bcm.Dividing by the volume in each cell gives the normalized reaction rates in tables 3 and 4. Note that other results in the literature [22,44] report the units as [reactions×10 −24 /source neutron].Figure 11 shows the 58 Ni reaction rates on a C/E basis.The (n,p) reaction rates show C/E values close to 1 particularly when using FENDL-3.2bdata, while the calculated (n,2n) reaction rates generally underpredict the experimental values regardless of which data library was used.These results are consistent with previous results in the literature [22,44], which also show a statistically significant underprediction of the (n,2n) reaction rates regardless of the choice of nuclear data library and MC code used for transport.
The SDR in the central cavity was calculated at each cooling time mentioned earlier ranging from 1 h to 60 d.At each cooling time, a source was defined for each activated cell with the energy distribution equal to the decay photon spectrum and the source strength equal to the decay photon intensity determined from the activated composition.The spatial distribution of each source is defined to be a uniform distribution over a bounding box of the cell, and cell rejection is then used to ensure that the spatial location of sampled photons is within the bounds of each activated cell.The total source strength, which by default serves as a multiplier on tally results in OpenMC, is equal to the sum of the decay photon intensities from each activated cell.The dose was calculated by running a photon transport calculation and tallying the flux in a spherical cell with radius 1.9 cm filled with air multiplied by  the ICRP-74 [53] flux-to-dose conversion factors 4 .Previous authors have shown that the choice of the dose conversion factors can have a major impact on the calculated dose [22]; for the present purposes, our goal is to show that the OpenMC results are in line with previous calculations and thus we have not performed an exhaustive comparison using dose conversion factors from different sources.As with the neutron transport simulations, photon transport simulations were run with 10 8 source particles divided into 100 batches.Simulations utilized ENDF/B-VIII.0photoatomic cross sections, which are based on EPICS2014 [54]. Figure 12 shows the SDR as calculated by OpenMC compared to the experimental measurements taken from the Geiger-Müller counter.Overall, the OpenMC calculated SDR agrees reasonably well with experimental values.A comparison of dose rates and C/E values at five specific cooling times, shown in table 5, puts our results in context with previous results.The calculated dose at short cooling times is generally within experimental uncertainties.At longer cooling times, there is an overprediction in the SDR relative to experimental values.However, this overprediction has been observed in previous calculations-for example, calculations made with MCR2S (MCNP and FISPACT) in 2010 [22].More recent results using MCR2S [44] with OpenMC and FISPACT-II had even higher C/E values, as seen in the bottom of figure 12, but this was attributed to the fact that the mesh used for activation overlapped with void in the cavity, resulting in a photon source

Conclusion
In this article we have demonstrated the validity of OpenMC for fusion-relevant, fixed-source transport simulations, activation calculations, photon source generation, and R2S SDR workflow management.The features developed herein to support SDR calculations represent an improvement in the state of the art for fusion neutronics and provide this community with the first fully open-source, nuclear analysis toolchain capable of SDR workflows.Consolidating the transport and activation calculations in a single software package significantly simplifies the SDR workflow and allows for more accurate solutions as the burnup matrix is generated directly from the continuous energy MC solution.
The cell-based R2S workflow that has been demonstrated here builds on OpenMC's existing capabilities while also relying on new features that have been added to the code specifically to support this workflow, including: decay photon source generation, cell rejection sampling, Walker alias sampling, enhancements to the IndependentOperator class, and support for decay photon spectra in depletion chain files.Each of the individual physics components in the R2S workflow has been verified through code-to-code comparisons.Neutron transport and photon transport simulations were performed in OpenMC and MCNP; the results show that, given the same nuclear data and problem description, the tallied flux spectra between the codes agrees within statistical uncertainty.Inventory calculations were carried out with OpenMC through the openmc.depletemodule and FISPACT-II.A comparison of predicted material compositions demonstrated excellent agreement, as seen in figure 5. Verification of the decay photon source generation was also performed by comparing to FISPACT-II.Although there are some subtle differences in how decay photon spectra are handled, for the common case of unstable nuclides with a discrete photon spectrum, the collapsed 24 energy group spectra calculated by OpenMC and FISPACT-II agree within a fraction of a percent (figure 7).
A limited validation of the full cell-based R2S workflow was carried out through a simulation of the first experimental campaign from the FNG ITER dose rate benchmark problem.For short cooling times, the dose calculated by OpenMC generally agrees with the experimental measurements accounting for uncertainties, shown in figure 12.For longer cooling times, an overprediction of the shutdown dose was observed relative to experiment, which is consistent with many previous studies.The reason for this overprediction is not known but could be attributed to many factors, including the choice of neutron cross sections, activation reactions and cross sections, decay data, photon transport data, and flux-to-dose conversion factors.The choice of the dose conversion factors by itself can result in a ±20% swing in the observed dose values (e.g.see the results in [22] wherein two different sets of dose coefficients were used).Sauvan et al [24] also made a crude estimate of the simulation uncertainties, accounting for material densities, cross sections, source power, isotope selection, irradiation time profile, and statistical errors and arrived at a total uncertainty of 17%.In conjunction with the 10% stated experimental uncertainty, it is clear that we should regard the simulation results as rough estimates and that it would be unreasonable to expect perfect agreement with experiment given the magnitude of the uncertainties.Further validation against other experiments, such as ones performed on JET [14], would provide further confidence in the code system and is recommended as future work.
While the scope of this work was limited to a cell-based R2S workflow for the sake of demonstrating the integration of physics solvers in OpenMC, we recognize that a cell-based workflow imposes constraints on the resolution of cells in order to obtain an accurate solution.In the case of the FNG dose benchmark problem studied here, sufficiently fine resolution in the number of cells was already present.However, for a general problem, a large spatial extent of individual cells may lead to poor resolution of the decay photon source and ultimately impact the accuracy of the solution.A mesh-based R2S workflow would alleviate the issues relating to spatial resolution and could naturally build upon the same software infrastructure that has been outlined here.A mesh-based workflow would also allow one to easily carry out a mesh resolution study to determine the necessary spatial fidelity required for a given problem, as has been demonstrated by others [55].
The end-to-end nuclear modeling and simulation capabilities developed and validated in OpenMC and described herein are well suited to both physics exploration for fusion technology development as well as for production-quality design iterations within fusion industry.Providing a single open-source software platform for both exploration and production analyses enables closer collaboration between researchers in academia and industry in a way that will facilitate the nuclear technology development necessary to design and build a functioning fusion power plant.Through a combination of accessibility, scalability, and simplicity, OpenMC and the capabilities developed in this work aim to enable fast, high-fidelity design iterations for fusion technology.

Figure 1 .
Figure 1.Key aspects of the R2S workflow that will be automated by OpenMC's Python API.Yellow rectangles represent physics modules and blue parallelograms represent intermediate data.

Figure 2 .
Figure 2. CAD model of the FNG dose benchmark showing the dimensions, materials, and cell segmentation for the constructive solid geometry (CSG) transport model.Highlighted in orange are two cells that are used in section 3 for the purpose of detailed comparison of neutron spectra and nuclide inventories between codes.One quarter of the model has been cut away for visualization purposes.

Figure 3 .
Figure 3. Multigroup flux tally comparison in cell 160 (a) and (b) and cell 373 (c) and (d) of the FNG dose benchmark between MCNP and OpenMC. )

Figure 4 .
Figure 4. Comparisons of the agreement between neutron transport solutions for OpenMC and MCNP for all 299 cells in the FNG dose benchmark.The histogram in (a) shows the probability distribution of relative differences between OpenMC and MCNP for all cell-energy bins pairs that have a converged solution according to MCNP (relative error <10%).The distribution shown in (b) displays the fraction of energy bins below a certain relative difference between OpenMC and MCNP for all cells.From this plot we glean that >80% of the energy bins in all cells show less than a 10% difference between OpenMC and MCNP and >95.7% have differences less than 20%.The bottom row depicts the cell-averaged relative error in each energy bin, showing excellent agreement between OpenMC and MCNP (c).The shaded region in (c) shows the bounds of the and maximum relative error for each energy bin across all cells in the model.The probability distribution shown in (d) demonstrates a high degree of similarity between the relative error distribution of all cell-energy bin pairs in both models.

Figure 5 .
Figure 5. Nuclide inventories for select nuclides for cell 160 (a) and cell 373 (b) are shown as calculated by OpenMC in coloredand by FISPACT-II in gray circular markers, showing excellent agreement for both long-and short-lived nuclides, as well as metastable states over a complex irradiation scenario.

Figure 6 .
Figure 6.Comparison of the total activity in all cells as computed by both OpenMC and FISPACT-II (a) as well as the relative difference between the two as a function of time (b).

Figure 7 .
Figure 7.Histogram of relative differences between OpenMC and FISPACT-II for the calculated average energy produced in each of the 24 energy groups from each unstable nuclide in ENDF/B-VIII.0 with a discrete spectrum.

Figure 8 .
Figure 8. Discrete photon spectra generated by OpenMC for an activated steel composition from the FNG dose benchmark (top) along with the same distribution binned onto a 24 energy group structure compared to FISPACT-II (bottom).

Figure 9 .
Figure 9. Dose per source photon in the central cavity of the FNG dose benchmark corresponding to a photon source in a single cell computed by OpenMC and MCNP (top) along with the relative difference between the two codes (bottom).Dashed lines correspond to twice the reported standard deviation from MCNP.
1 h, 6 h, 12 h, 16 h, 20 h, 1 d, 2 d, 3 d, 4 d, 5 d, 7 d, 9 d, 12 d, 15 d, 18 d, 21 d, 30 d, and 60 d.While computational results are typically compared at cooling times of 1 d, 7 d, 15 d, 30 d and 60 d, the expanded list of cooling times allows us to assess whether the dose rate at short cooling times matches the experimental values as well.

Figure 10 .
Figure 10.Table of nuclides showing which nuclides are included in the full and reduced depletion chains based on ENDF/B-VIII.0.

Figure 11 .
Figure 11.Ratio of calculated (C) to experimental (E) reaction rates in 58 Ni in the FNG dose benchmark.

Figure 12 .
Figure 12.Campaign 1 experimental and OpenMC-calculated shutdown dose rate in the central cavity in [Sv/h] (top) and ratio of calculated (C) to experimental (E) shutdown dose rate (bottom).

Table 1 .
Summary of SDR workflows from the literature a .

Table 2 .
Irradiation schedule for the FNG shutdown dose rate experiment, campaign 1.

Table 3 .
Experimentally measured and calculated 58 Ni(n,p) 58 Co reaction rate in [reactions/cm 3 -source] for a hypothetical 1 atom/b-cm concentration of 58 Ni.