Dear Paul:
I am trying to produce a plot of the normalized neutron flux (per unit lethargy) as a function of neutron energy for the reactor described in the attached ,pynb file. Also attached is a generic .jpg file of what I’m expecting.
I have trouble working with arrays in python. Can you help me with this? Is there a simple script for such a plot?
Sincerely,
Shawn Wachter
MCFR-RPU-final-001.ipynb (37.3 KB)
Hello, I am trying to generate a neutron flux plot. Can you please guide me how you did the plot for neutron flux. I am trying to generate flux tally with the below code
<?xml version="1.0"?>regular
50 1
<lower_left>-70.71 -70.71</lower_left>
70.71 70.71
1
1
flux
where 70.71 is the radius of my reactor. Now, In the tally.out file I got something
============================> TALLY 1 <============================
Mesh Index (1, 1)
Total Material
Flux 12.2047 +/- 6.30332E-02
Mesh Index (2, 1)
Total Material
Flux 11.7069 +/- 6.12077E-02
…
for the first two mash. And the rest is 0. I want to plot the neutron flux from center to radius of the reactor. I am not getting such values. Am I missing something? If I had the value I would use them to generate a downward plot in Microsoft excel. But, so far I am stuck on not getting the proper data in flux tally. Please help.
Hi Shawn,
My apologies for the slow response on this. I’ve put together a quick example notebook showing how one would tally a flux spectrum using the Python package. Let me know if you have any questions about it.
Best regards,
Paul
Hi Sharif,
The problem in your example is that the specified for the mesh is actually the width of a single mesh element. This means that the width of the your entire mesh is actually 50*70.71. Instead you’ll probably want to specify the width as 1.4142 70.71.
Best regards,
Paul
Hello, Dr Paul. Thank you for your reply. I understand my problem. Now, I get the deserved result in the tally.out. Is there any way to plot the data in tally.out with openmc?
Thanks and regards,
Sharif Abu Darda, Mohammad Talhi.
50 1 mash.rtf (6 KB)
When you run OpenMC, you also get a “statepoint” file which uses the binary HDF5 format. The statepoint file holds all the tally results that you find in tallies.out. To get tally results from the statepoint file and plot them, we generally recommend writing a Python script because 1) OpenMC includes its own Python package that allows you to interact with the various input/output files and 2) there are many popular third-party packages for plotting. To get an idea of what you can do with the Python package, I would recommend having a look at our example notebooks.
Best regards,
Paul
Hello, Dr, I have never used Python API before. I am on the Mac OS system. I open the Python API from the miniconda3/python/Contents/MacOS/python. Now, I believe this points to my root directory. Now If I put the “statepoint.50.h5” of my project in there and, from your guide, I can load the file in python API "
sp = openmc.StatePoint('statepoint.50.h5')", Can you please guide me on the commands to run on the pythonAPI after opening it to get the flux plot? I will load the statepoint file manually by above process.
A list of commands for this perticular query would be much helpful. I am tying to learn to work with python API also.
Thanks you.
Regards,
Sharif Abu Darda, Mohammad Talhi
Hi Sharif,
Once you have a statepoint file opened, it has a ‘tallies’ attribute that gives you access to tally metadata and results. Since in your problem you have a single tally (ID=1), to get that tally through the Python API, you could run:
sp = openmc.StatePoint(‘statepoint.50.h5’)
tally = sp.tallies[1]
flux = tally.mean.ravel()
Running tally.mean gives you a multidimensional numpy array where the dimensions correspond to filter combinations and scores. Since you only have a single filter (mesh) and a single score (flux), it’s easier to deal with a one-dimensional array, which you can get by running the ravel() method. To plot the flux, you could use the matplotlib Python package with something like:
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-70.71, 70.71, 50)
plt.plot(x, flux)
plt.xlabel(‘x position [cm]’)
plt.ylabel(‘flux’)
plt.show()
Hope this helps!
Best,
Paul
Dear Paul,
Continuing with this post, I have some questions.
I am simulating a single fuel element and, in my report, I would like to include a graph showing the flux spectrum. To get that spectrum, I created the flux tally following the steps of this example notebook.
Hi Javier,
The reason that log10 is used when calling logspace just has to do with how logspace treats its arguments. The bounds are treated as powers of ten so logspace(1, 3) would give you points between 10 and 1000. Hence, if you want points between E1 and E2, you need to pass log10(E1) and log10(E2) so that the exponential and log cancel out.
Otherwise, the normalization that you show looks correct to me.
Best,
Paul
Thanks Paul for your answer.