Best way to create DPA spatial distribution map


I want to get a spatial DPA (which can be computed from energy-binned flux tallies) distribution. Ideally, I could get this in the format of (x, y, z, flux) or (x, y, z, DPA) per element.

What is the best way to do this?

I’m currently anticipating @pshriwise’s PR which will enable the generation of a .vtk with more than just one filter as well as write elemental coordinates to the statepoint file.

Although it will be written to and easily visualizable with the .vtk, I suspect I still need to extract the data from the statepoint file to manipulate the data.

When following the method as described in this tutorial, I found that with a mesh of about a million elements, just binning by two energy filters took python about 3 hours to complete, specifically these lines of code:

    thermal_flux = tally.get_values(scores=['flux'], 
                                    filter_bins=[((0.0, 1.e+06),)])
    fast_flux = tally.get_values(scores=['flux'],
                                 filter_bins=[((1.e+06, 5.e+06),)])

The data that I have for DPA for graphite is binned into 100 energies. Since the above code parses through the statepoint sequentially, to run this binning over 100 energies would take 50*3 = 150 hours.

Is there a better way to do this?

Well, 3 hours for a single call to get_values doesn’t sound good! It certainly seems like something is not scaling well within that method, but after looking at the code, nothing is really jumping out at me as being the obvious culprit. Would you be able to run a quick profile on a call to that method? All you would need to do is:

  1. Create a few-line .py script that loads a statepoint and then makes a single call to tally.get_values as above.
  2. Call the script with python -m cProfile
  3. Paste the results here

That would be really helpful to narrow down what’s going on.

Also, for calculating DPA, you may want to specify the damage-energy score. This is based on the MT=444 data that is computed by NJOY/HEATR and is the typical method that I’ve seen people use to calculate DPA from Monte Carlo simulations.

If there is a tally.out file,just use regular expression to read and filter data.Recently i try to read a complex geometry reactor flux tally statepoint.h5 file,I also failed cause memory overflows,and it will cause too much time.

I ran a profile on the single call to tally.get_values for my mesh, and it took 4 hours and the backtrace was so long it didn’t fit in the terminal. I should’ve piped it to a file, so now it’s running again… update in a few hours.

@paulromano I ran the profile again, and it was ~5000 lines. I saved it as an .xml so I could upload it here.

Also, do you have any references for the damage-energy method you mention? I’d like to look into that.

profile.xml (396.1 KB)

Thanks for the profile @kellythomas! I’ll have a look and see if I can figure out what’s going on.

Regarding damage-energy, a good place to start would be the NJOY user’s manual, in particular the section on the HEATR module (which is responsible for generating this data). I also gave a presentation recently at Argonne on simulating radiation damage, which I think you’ll find helpful. You can find the slides here.

Ok, so looks like the super slow get_values call is due to an issue that has already been fixed by @pshriwise but just hasn’t hit a release yet:

In our next release of OpenMC, the calls you are making to get_values should take a trivial amount of time.

1 Like

Thank you @paulromano! Do you know when the next release of OpenMC might be?

I’m hoping to get a release out in the next week or two. We’re just wrapping up a few final issues. Stay posted!