plot of normalized neutron flux per unit lethargy as a function of neutron energy

[shawn.wa...@gmail.com]

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

[Sharif Abu Darda]

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"?>

<tallies>

    

    <mesh id="1">

        <type>regular</type>

        <dimension>50 1</dimension>

        <lower_left>-70.71 -70.71</lower_left>

        <width>70.71 70.71</width>

    </mesh>

    

    <filter id="1" type="mesh">

        <bins>1</bins>

    </filter>

    

    <tally id="1">

        <filters>1</filters>

        <scores>flux</scores>

    </tally>

    

</tallies>

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.

[Paul Romano]

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

--
You received this message because you are subscribed to the Google Groups "OpenMC Users Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to openmc-users+unsubscribe@googlegroups.com.
To post to this group, send email to openmc...@googlegroups.com.
Visit this group at https://groups.google.com/group/openmc-users.
For more options, visit https://groups.google.com/d/optout.

[Paul Romano]

Hi Sharif,

The problem in your example is that the <width> 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 <width>1.4142 70.71</width>.

Best regards,

Paul

--

[Sharif Abu Darda]

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.

[Paul Romano]

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

[Sharif Abu Darda]

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

[Paul Romano]

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

To unsubscribe from this group and stop receiving emails from it, send an email to openmc-users...@googlegroups.com.

[Javier Gonzalez]

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.

--------------------------------

# Create equal-lethargy energies to put in filter
energies = np.logspace(np.log10(1e-5), np.log10(20.0e6), 501)
energy_filter = openmc.EnergyFilter(energies)

# Create tally with energy filter
tally = openmc.Tally(tally_id=100, name=’Flux’)
tally.scores = ['flux']
tally.filters = [energy_filter]

--------------------------------

Now, I want to transform the y-axis units. Specifically, in the y-axis, I want the units “flux per unit lethargy [neutrons/cm²s]” and, in the x-axis, “Energy [eV]”. My questions are:

1- Why the equal-lethargy energies to put in the filter were created considering “np.log10(x)”? Isn’t the lethargy calculated using the “np.log(x)” function?

2- To transform the flux units from “flux [neutrons-cm/source neutron]” to “flux per unit lethargy neutrons/cm²s]”, should I multiply the tally by the normalization factor (presented in another post) and divided it by the lethargy? Something like this:

--------------------------------

t = sp.get_tally(name='Flux')

flux = t.mean.ravel()

factor = (P*nu_fission)/(fission_rate*vol*Q*k_eff)

delta_u = np.log(energies[1:501]) - np.log(energies[0:500])

flux1 = flux*factor/delta_u

plt.semilogx(energies[:-1], flux1, color=”red”)

plt.xlabel(“Energy [eV]”)

plt.ylabel(“Flux per unit lethargy [neutrons/cm$^2$s]”)

--------------------------------

Thank you,

Javier

On Friday, April 27, 2018 at 7:35:42 AM UTC-4, Paul Romano wrote:

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

On Tue, Apr 10, 2018 at 10:18 AM, <shawn.w...@gmail.com> wrote:

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

--
You received this message because you are subscribed to the Google Groups "OpenMC Users Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email to openmc...@googlegroups.com.

[Paul Romano]

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

To unsubscribe from this group and stop receiving emails from it, send an email to openmc-users...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/openmc-users/2f14f6ca-2470-4f75-972a-977e9c444e42%40googlegroups.com.

[Javier Gonzalez]

Thanks Paul for your answer.

To unsubscribe from this group and stop receiving emails from it, send an email to openmc...@googlegroups.com.