I am running a simple activation analysis comparison (for Fixed neutron Source Transmutation Calculations) case by benchmarking openmc results against FISPAC-II code to confirm that the two codes give identical results.
The case that I am running is a sphere of pure Manganese (Mn) with 14MeV source at the centre of the sphere. Please see attached input file for the run for openmc Mangane_deplet.py (2.1 KB)
.
My mental model is that Manganese has only one stable nuclide, Mn-55, with an abundance of 100%. Hence, Mn-56 should be the main activation product if I irradiate it with a flux of 2.3760E+15 neutrons/cm2/s (i.e., cm-2 s-1) for 6.307e+7s and allowed to decay for only 1 seconds. Since the half-life of Mn-56 is about 2.578 hours.
The FISPAC-II result is inline with my expectation and only Mn-56 isotope is the main activity. However, the openmc results is not inline with my expectations but rather it produced Mn-54 as the only main activity. I have checked the openmc nuclear data library and I found that the nuclear data library for Mn-56 nuclide is missing. I then checked different nuclear data libraries for openmc (both h5 and ACE files) and in all cases the data library for Mn-56 nuclide is missing.
My question is why is openmc getting it wrong and why is the Mn-56 nuclide library missing as this data exist for the FISPAC-II code. Kindly have a look at my attached input and let me know if I have got the input specification wrong for openmc.
Thank you so much for your help.
Hi @Daniel thanks for bringing this up it’s super important! There are a number of issues going on here:
You’re right that Mn56 isn’t present in the main incident neutron libraries, which is strange, but I think is the main cause of the problem perhaps?.
I’m not sure which depletion chain you’re using, but the endf ones on openmc.org have Mn56 in them.
FISPACT doesn’t give you any real choices about the time integration scheme, whereas OpenMC gives you many options. Under the hood, FISPACT essentially inserts a bunch of intermediate timesteps to limit the error due to the time integration scheme. Right now, OpenMC doesn’t do that. If you use higher order integrators you should get better agreement with FISPACT, right now you’re just using the PredictorIntegrator.
@Daniel In your script I don’t actually see any logic for getting the activity. When I add the following code to get the activity by nuclide (using an ENDF/B-VII.1 depletion chain):
I believe the main issue was to do with the way the results were being postprocessed because I did use Paul’s above python code to postprocess the results and had similar results as quoted in his post despite missing the Mn56 neutron library data.
I did used most of the various depletion chain found at Depletion Chains | OpenMC with my postprocessing script but it did not resolved the issue until I used Paul’s above python code to postprocess the results.
The results from the two codes (FISPACT and OPEN) now seems to be heading towards the same direction of travel when Paul’s above python code is used to postprocess the results. My next step is to proceed to use the higher order integrators within Openmc to achieve even more closer agreement between the to codes. I will keep you updated with my findings.
Paul, with regards to your queries on how the results was postprocessed, I did used the depletion plotter script from the Silver example written by Shimwell (plot_activity_vs_time.py (817 Bytes)) and also a slightly modified version to suite my own need (Myplot_activity_vs_time.py (2.6 KB)). Note that I was using the latest openmc-depletion-plotter 0.3.1 with the latest openmc 13.3. I will proceed with the rest of my benchmark analysis and keep you posted.
Once again thank you both for your help.
An important thing that is going on here is that @Shimwell’s plot_activity_vs_time calls the export_to_materials method of the openmc.deplete.Results class which means that nuclides that don’t have transport cross sections (i.e. Mn56) won’t be written out and therefore won’t be plotted. If you want to properly get the activity of the material, use Paul’s suggestion above. The activated_mat variable he uses will have all the nuclides produced during the depletion calculation.
Nice spot Ethan. I didn’t realize export_to_materials was excluding nuclides that are not present in the cross_sections.xml. I see it in the doc string so I must have just missed it first time .
Paul’s method does look better.
I’ve updated the openmc-depletion-plotter package and version 0.3.2 should fix this issue.
pip install openmc_depletion_plotter==0.3.2
I guess @Daniel was using the graphical user interface with the package, hence no actual extracting of nuclides in the provided script?
@Shimwell thank you for fixing the issue in the openmc_depletion_plotter package. I have now installed the 0.3.2 version and can confirm that it is working as expected. I was using the graphical user interface with the package and also in another instant I had to hack your code in order to be able to list all the nuclides being tracked in a tabulated format.
@eepeterson and @paulromano as promised below is the summary of result of following on analysis based on your above expert advice:
@eepeterson as you did recommended, I have extended my analysis to test each of the advanced depletion integrators implemented in the OpenMC code in turns i.e. predictor, CE/CM(CECMIntegrator), CE/LI(CELIIntegrator), CF4(CF4Integrator), EPC-RK4(EPCRK4Integrator), LE/QI(LEQIIntegrator), SI-CE/LI(SICELIIntegrator) and SI-LE/QI(SILEQIIntegrator). All the other integrators, except for predictor, provided similar results, hence, in the below summary comparisons of results with FISPACT I will only dwell on the CE/CM integrator results to make comparisons with the FISPACT code.
Summary of Analysis
The setup consisted of a sphere of radius 100cm with 14MeV neutron source at its centre. The sphere was irradiated for 6.30700E+07 secs using a flux of 2.37600E+15 n/cm**2/s. Followed by a cooling time of 45000 secs (i.e. 5 half-lives of Mn56). Paul’s (@paulromano) script was modified slightly (by setting last_step = results[-2] to extract result for Step 1 and last_step = results[-1] to extract result for Step 2) so that I can extract the result for the two steps individually. The OpenMC runs used TENDL-2019 nuclear data library whereas the FISPACT runs used TENDL-2017 nuclear data library (see script used for the run Mangane_deplet.py (2.9 KB)). The results were obtained as follows:
Step 1 : After irradiating for time = 6.3070E+07 SECS and flux =2.3760E+15 /cm^2/s
For the OpenMC runs I was missing some nuclear libraries including the Mn56 neutron library data, hence, it was unable to have the same pathway as that of the FISPACT case which had all the nuclear data libraries present. The pathway reported by FISPACT was: Mn 55 ---(R)--- Mn 56 ---(d)--- Fe 56 ---(R)--- Fe 57 ---(R)--- Fe 58 ---(R)--- Fe 59 ---(d)--- Co 59 ---(R)--- Co 60 ---(L)---.
The dominant nuclides listed by the OpenMC code included Mn54 nuclide which I believe results from its (n, 2n) reaction threshold whereas the FISPACT code did not list the Mn54 nuclide which I am a bit puzzled. Any taught on these findings will be much appreciated and whether there is a future plan to release a more detailed nuclear libraries for the OpenMC code which will improve the accuracy of the results between FISPACT code.
Thank you.
@Shimwell Thank you for your comments. As suggested, I have emailed my results to one of the FISPACT experts developer/user at UKAEA. Also, I will raise it on their forum as well. I do agree that using exactly the same TENDL nuclear data libraries for both codes will provide consistency in calculation, which I will do in my next analysis. Nonetheless, I do not expect it to improve the results significantly as I believe the nuclear data libraries do not change significantly between versions. I will keep you updated once I have completed my analysis.
Hello @eepeterson, @paulromano and @Shimwell,
I have had the opportunity to enhance my modelling with your support, as well as the assistance from FISPACT experts and developers/users at UKAEA. Currently, both codes report very similar nuclides, with the exception that OpenMC activity is approximately two orders of magnitude lower than that of FISPACT. Please refer to the attached image for your perusal.
The FISPACT experts and developers/users at UKAEA have confirmed that their results are correct and broadly align with rough hand analytical calculations. Consequently, I am uncertain why the OpenMC activity results are lower.
Both codes utilised equivalent nuclear data, and the simulations commenced with the same number of atoms (Mn-55 = 1.09617E+25 atoms) before undergoing irradiation for the same duration. FISPACT-II, being an inventory code, lacks spatial dimensionality and operates solely in the dimension of time. Since both codes resulted in the same number of atoms before irradiation, it gives me confidence that the geometry setup in OpenMC is of a similar volume to that of FISPACT-II.
I attempted using both the Discrete 14e+6eV and Tabular source with flux probability energy groups from FISPACT-II in OpenMC, and both resulted in identical activity, which remains lower than that of the FISPACT-II activity results.
I would appreciate it if you could kindly review my input file again Mangane_deplet.py (2.6 KB)
to confirm if there is any reason why the OpenMC activity results would be lower than that of the FISPACT case.
Thank you all for your continuous help and support. This is really much appreciated.
Hi Daniel welcome back glad to see you are working on this, there has been some development on related work in the last 7 months
@nuclearbae has been working on making a transport free and therefore more similar to fispact method in OpenMC. It is not finished by you can see the Pull request here and an example here (Jin has made some good comments on that issue that you might be interested in).
I should mention there is an open source inventory code ALARA that could also be compared and also ORIGEN if you want some more checks.
For your example I would guess that because OpenMC is transporting the neutrons through a geometry where you have a point source in the centre there are going to be lots that escape the small 3cm radius sphere and leak out the sides. I wonder if the source or geometry can be changed to match the approximations fispact is using. Would we get the same difference if the sphere wall 300cm radius for example?
@Daniel The source rate you specify for the integrator object is the total neutron source rate for the problem and not the flux. For your problem, OpenMC is doing transport in a sphere so the resulting flux that you will get will be lower than the 2.3760E15 source rate you specified (by roughly a factor of 100 I suspect). If you tally the flux spectrum in your sphere with the same energy bin structure that you use to supply the flux for FISPACT then that will tell you what your source rate needs to be to get matching results. It would be helpful to post your FISPACT inputs as well if you have them.
Ah yes that makes sense, thanks Ethan. The flux input to fispact is not the same thing as the energy distribution of the point source.
The flux used in fispact was properly a flux spectra tally made during a transport simulation but these get seperated and these flux files are distributed without the original geometry, materials, settings that generated them.
So if you want a direct equivalent I think Ethans suggests is the way to go. Add a flux spectra tally to your OpenMC transport and that will give you the flux sprctra to input into inventory codes like fispact.
Reverse engineering the original source term definition that results in the flux tally you are currently inputting into fispact might not be possible without knowing more about how it was obtained (geometry, materials, source, settings etc).
@Shimwell and @eepeterson, many thanks for your swift responds. I will adjust my approach accordingly. I will repeat my simulations but this time I will collect the neutron tally in OpenMC and pass to FISPACT-II as suggested. This was also suggested by the FISPACT experts and developers/users at UKAEA but I am yet to perform this run.
@eepeterson as requested please find attached the FISPACT input file used for the run for your perusal Manganese_1Kg.i.pdf (111.1 KB)
I have converted this into pdf in order to be able to attached to this forum.
@eepeterson and @Shimwell, I have tally the flux spectrum in the sphere with the same energy bin as the flux energy group structure in FISPACT and extracted the results. However, I am unsure whether to use the actual flux spectra tally (i.e. the mean flux score) values as it is which will be in unit of particle-cm/source and pass it on to FISPACT or I have to convert the tally scores into a more natural flux unit like particles/cm2-s by multiplying by the source rate which I don’t know. Thank you.
@Shimwell and @eepeterson I appreciate your assistance and support. Following the adjustments to my model as per your recommendations, I have achieved satisfactory agreement between OpenMC and FISPACT-II, particularly noticeable in the results for Mn as illustrated below.
However, I encountered a significant discrepancy in a specific case involving Na material as shown below:
The differences between OpenMC and FISPACT-II results were considerably larger. Upon further investigation, I discovered that OpenMC Version 0.13.3 did not accurately model Na. For instance, when I set the flux to 0.0 n/cm^2/s, assuming no neutron irradiation, and simulated the run for the same duration, I expected to end the simulation with the same number of atoms (2.61948E+25) of Na23_stable nuclide that I started with. To my surprise, the simulation resulted in Mg23 with an activity of 1.24715E+01 (Bq) for the initial 6.3070E+07 seconds, and for the subsequent cooling time of 45000 secs, it resulted in Mg23 with an activity of 3.75345E+05 (Bq).
I would appreciate it if you could kindly explain why OpenMC’s behaviour deviated from my expectations. I have attached my OpenMC input script and its output below for your review. Sodium_deplet.py (2.5 KB) Output_Sodium.pdf (110.9 KB)
I would be curious to know if or how you considered self-shielding in FISPACT (keywords: PROBTABLE or SSFGEOMETRY)?
My experience so far is that when using a properly tallied spectrum in the 709 or 1102 group from OpenMC in FISPACT as an input for the activation of Zircaloy-4 (FLUX mode) that the best results are those not applying self-shielding corrections in FP.
Using the energy-based self-shielding correction (PROBTABLE) on top of the tallied neutron spectrum with dips results in a strong underprediciton of activity in FP (over-shielding).
Also, regarding the TENDL libraries I noted differences between the years. See below for example a final result comparison in FP for TENDL2019 vs TENDL2017 at 600K.
@KUNG Thank you for your detailed response and insights. I appreciate you sharing your observations.
For my modelling in FISPACT, I did not apply any energy-based self-shielding corrections, such as PROBTABLE or SSFGEOMETRY. It’s valuable to hear your experience with tallied spectra and the effects of self-shielding corrections.
Regarding the differences between TENDL2019 and TENDL2017, thank you for highlighting these. It’s an area I’ll need to explore further.
Currently, I am in the process of performing experimental measurements on a known material. The plan is to benchmark these experimental results against predictions from both FISPACT and OpenMC. This should provide clarity on which approach or configuration yields the most reasonable predictions.
I will share any findings that might contribute to this discussion once the benchmark is complete.
.
