adding inventory code example showing deplete/activate/transmute with pre-calculated flux spectrum

[shimwell]

It would be great to add an example of irradiation with a known/pre-calculated flux and deplete/activate/transmute a material.

This example would be similar to inventory codes (ALARA, ORIGEN and FISPACT) and would not perform neutron transport.

While OpenMC can do more advanced deplete/activate/transmute with neutron transport these inventory codes don't account for transport.

So if we want a simplified deplete/activate/transmute example to compare OpenMC with inventory codes we can "downgrade" the deplete/activate/transmute methods available in OpenMC to make this simplified inventory simulation.

[shimwell]

example python code, currently needs this branch
https://github.com/jbae11/openmc/tree/microxs_from_mg_flux
which might be merged in via this openmc PR openmc-dev/openmc#2755

# This script simulates depletion of a material with a known neutron flux spectrum
# This is the OpenMC equivilent to inventory codes like ALARA, ORIGEN and FISPACT

import openmc
import openmc.deplete
import math

# users might want to change these to use specific xml files to use particular decay data or transport cross sections
# the chain file was downloaded with
# pip install openmc_data
# download_endf_chain -r b8.0
# openmc.config['chain_file'] = 'chain-endf-b8.0.xml'
# openmc.config['cross_sections'] = 'cross_sections.xml'

# makes a simple material from Silver
my_material = openmc.Material() 
my_material.add_element('Ag', 1, percent_type='ao')
my_material.set_density('g/cm3', 10.49)
my_material.depletable=True

# volume must set the volume as well as openmc calculates number of atoms
my_material.volume = 100 # units cm3
materials = openmc.Materials([my_material])
# We define timesteps together with the source rate to make it clearer
# units are seconds for time and number of neutrons for source rates
timesteps_and_source_rates = [
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),  # should saturate Ag110 here as it has been irradiated for over 5 halflives
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
]

timesteps = [item[0] for item in timesteps_and_source_rates]
source_rates = [item[1] for item in timesteps_and_source_rates]

flux_in_each_group=[1.1e-7, 1.2e-6, 1.3e-5, 1.4e-4]

micro_xs = openmc.deplete.MicroXS.from_multi_group_flux(
    energies='CCFE-709',
    multi_group_flux=[0.1]*709,
    temperature=293,
    chain_file= openmc.config['chain_file']
)

# constructing the operator, note we pass in the flux and micro xs
operator = openmc.deplete.IndependentOperator(
    materials=materials,
    fluxes=flux_in_each_group,
    micros=[micro_xs],
    reduce_chain=True,  # reduced to only the isotopes present in depletable materials and their possible progeny
    reduce_chain_level=5,
    normalization_mode="source-rate"
)

integrator = openmc.deplete.PredictorIntegrator(
    operator=operator,
    timesteps=timesteps,
    source_rates=source_rates,
    timestep_units='s'
)

# this runs the depletion calculations for the timesteps
# this does the neutron activation simulations and produces a depletion_results.h5 file
integrator.integrate()

results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")

times, number_of_Ag110_atoms = results.get_atoms(my_material, 'Ag110')

for time, num in zip(times, number_of_Ag110_atoms):
    print(f" Time {time}s. Number of Ag110 atoms {num}")

import matplotlib.pyplot as plt
plt.plot(times, number_of_Ag110_atoms)
plt.xlabel('Time [s]')
plt.ylabel('number of Ag110 atoms')
plt.show()

[shimwell]

I have updated the example to make the micro_xs from nuclides that are in the materials instead of the whole chain file. This is much quicker as in this example only 2 nuclides are present and the chain file contains hundreds.

# This script simulates depletion of a material with a known neutron flux spectrum
# This is the OpenMC equivilent to inventory codes like ALARA, ORIGEN and FISPACT

import openmc
import openmc.deplete
import math

# users might want to change these to use specific xml files to use particular decay data or transport cross sections
# the chain file was downloaded with
# pip install openmc_data
# download_endf_chain -r b8.0
# openmc.config['chain_file'] = 'chain-endf-b8.0.xml'
openmc.config['chain_file'] = '/home/jshimwell/ENDF-B-VIII.0-NNDC/chain-nndc-b8.0.xml'
# openmc.config['cross_sections'] = 'cross_sections.xml'
openmc.config['cross_sections'] = '/home/jshimwell/ENDF-B-VIII.0-NNDC/h5_files/cross_sections.xml'

# makes a simple material from Silver
my_material = openmc.Material() 
my_material.add_element('Ag', 1, percent_type='ao')
my_material.set_density('g/cm3', 10.49)
my_material.depletable=True

# volume must set the volume as well as openmc calculates number of atoms
my_material.volume = 100 # units cm3
materials = openmc.Materials([my_material])
# We define timesteps together with the source rate to make it clearer
# units are seconds for time and number of neutrons for source rates
timesteps_and_source_rates = [
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),  # should saturate Ag110 here as it has been irradiated for over 5 halflives
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
]

timesteps = [item[0] for item in timesteps_and_source_rates]
source_rates = [item[1] for item in timesteps_and_source_rates]

flux_in_each_group=[1.1e-7, 1.2e-6, 1.3e-5, 1.4e-4]

micro_xs = openmc.deplete.MicroXS.from_multi_group_flux(
    energies='CCFE-709',
    multi_group_flux=[0.1]*709,
    temperature='294', # endf 8.0 has ['1200K', '2500K', '250K', '294K', '600K', '900K']
    chain_file= openmc.config['chain_file'],
    nuclides=my_material.get_nuclides()
)

# constructing the operator, note we pass in the flux and micro xs
operator = openmc.deplete.IndependentOperator(
    materials=materials,
    fluxes=flux_in_each_group,
    micros=[micro_xs],
    reduce_chain=True,  # reduced to only the isotopes present in depletable materials and their possible progeny
    reduce_chain_level=5,
    normalization_mode="source-rate"
)

integrator = openmc.deplete.PredictorIntegrator(
    operator=operator,
    timesteps=timesteps,
    source_rates=source_rates,
    timestep_units='s'
)

# this runs the depletion calculations for the timesteps
# this does the neutron activation simulations and produces a depletion_results.h5 file
integrator.integrate()

results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")

times, number_of_Ag110_atoms = results.get_atoms(my_material, 'Ag110')

for time, num in zip(times, number_of_Ag110_atoms):
    print(f" Time {time}s. Number of Ag110 atoms {num}")

import matplotlib.pyplot as plt
plt.plot(times, number_of_Ag110_atoms)
plt.xlabel('Time [s]')
plt.ylabel('number of Ag110 atoms')
plt.show()

[jbae11]

I might be wrong (and in a lot of cases am), but wouldn't doing nuclides=my_material.get_nuclides() have the risk of not taking into account secondary nuclides (fission / activation products, decay daughters)? I'll probably test this and see if it makes a difference for various cases, but looks great! I'll try this

[shimwell]

Yep I think you are totally right about that. Only nuclides and their direct reaction products and their decay products are tracked with this approach. It would be nice to have something that could get all the surrounding nuclides to an arbitrary thickness.

Something like
nuclide + p
Nuclide + n
Nuclide + n+p
Nuclide - p
Etc

Shall I make a small function

[jbae11]

hmm I can imagine something like a recursive function to get all the chains and resulting nuclides from a 'fresh composition' - I think it'll be a good function in the Chain class right?

[shimwell]

I like the sound of that, perhaps a separate PR for a new method on the chain class

[shimwell]

Example corrected, thanks @jbae11 I had misunderstood flux arg in the IndependentOperator

# This script simulates depletion of a material with a known neutron flux spectrum
# This is the OpenMC equivilent to inventory codes like ALARA, ORIGEN and FISPACT

import openmc
import openmc.deplete
import math

# users might want to change these to use specific xml files to use particular decay data or transport cross sections
# the chain file was downloaded with
# pip install openmc_data
# download_endf_chain -r b8.0
openmc.config['chain_file'] = '/home/j/chain-endf-b8.0.xml'
# openmc.config['cross_sections'] = 'cross_sections.xml'


# makes a simple material from Silver
my_material = openmc.Material() 
my_material.add_element('Ag', 1, percent_type='ao')
my_material.set_density('g/cm3', 10.49)
my_material.depletable=True

# volume must set the volume as well as openmc calculates number of atoms
my_material.volume = 100 # units cm3
materials = openmc.Materials([my_material])
# We define timesteps together with the source rate to make it clearer
# units are seconds for time and number of neutrons for source rates
timesteps_and_source_rates = [
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),  # should saturate Ag110 here as it has been irradiated for over 5 halflives
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 1e20),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
    (24, 0),
]

timesteps = [item[0] for item in timesteps_and_source_rates]
source_rates = [item[1] for item in timesteps_and_source_rates]

# this is the precalculated neutron spectra in units of n/cm2
flux_in_each_group=[0.1]*709

micro_xs = openmc.deplete.MicroXS.from_multigroup_flux(
    energies='CCFE-709',
    multigroup_flux=flux_in_each_group,
    temperature=294, # endf 8.0 has ['1200K', '2500K', '250K', '294K', '600K', '900K']
    chain_file= openmc.config['chain_file'],
    nuclides=my_material.get_nuclides()
)

# constructing the operator, note we pass in the flux and micro xs
operator = openmc.deplete.IndependentOperator(
    materials=materials,
    fluxes=[sum(flux_in_each_group)*my_material.volume],  # Flux in each group in [n-cm/src] for each domain
    micros=[micro_xs],
    reduce_chain=True,  # reduced to only the isotopes present in depletable materials and their possible progeny
    reduce_chain_level=5,
    normalization_mode="source-rate"
)

integrator = openmc.deplete.PredictorIntegrator(
    operator=operator,
    timesteps=timesteps,
    source_rates=source_rates,
    timestep_units='s'
)

# this runs the depletion calculations for the timesteps
# this does the neutron activation simulations and produces a depletion_results.h5 file
integrator.integrate()

results = openmc.deplete.ResultsList.from_hdf5("depletion_results.h5")

times, number_of_Ag110_atoms = results.get_atoms(my_material, 'Ag110')

for time, num in zip(times, number_of_Ag110_atoms):
    print(f" Time {time}s. Number of Ag110 atoms {num}")

import matplotlib.pyplot as plt
plt.plot(times, number_of_Ag110_atoms)
plt.xlabel('Time [s]')
plt.ylabel('number of Ag110 atoms')
plt.show()

• edit updated api method name from_multi_group_flux to from_multigroup_flux

[jbae11]

okay so I've been comparing the results with some ORIGEN results (that match with the FISPACT results from the FNS benchmark) and I found that:

• we have to normalize the flux array that goes into the MicroXS.from_multi_group_flux function - that way the microxs values match those calculated by ORIGEN
• The units check out on the IndependentOperator flux argument, but I'm not sure if that has to normalized as well. I feel like the scaling is weird, but I'll keep testing to see what yields the most similar results (it definitely has a huge impact on the result).
• Sometimes the results are unrealistic (total number of atoms increases), so maybe there's something fundamental going on.

Given all this, do you know if we have a very well-documented benchmark we can follow for this workflow? like:

• known MG cross section
• known material with volume
• known collapsed XS
• known nuclide composition over time, etc.

[jbae11]

I've been playing around with the FNS data and this is what I have so far, it's obviously very wrong, but I'm not sure if it is due to some scaling error, or instabilities in the solver.. Any guesses? I checked the microxs values and the reaction yields and it looks to be similar to that of ORIGEN, so my guess is that the errors are from the oeprator and integrator side

[jbae11]

implementation here

[shimwell]

Is there any chance of adding time steps to the openmc one. Extra time steps between the existing time steps. I think something fispact does that openmc doesn't do yet is add additional internal timesteps if the user timesteps are too large. I think max timestep can be found from the stiffness of the matrix automatically but for now perhaps add a few more user tomesteps to the openmc one

Perhaps add extra args to the IndependentOperator

reduce_chain=True,
reduce_chain_level=5

I suspect that if you plot the numbers of different nuclides then there are some spikes for individual nuclides

[jbae11]

That was a great idea - looks like a lot of the craziness disappeared. I'm probably still messing up something with scaling cause it still looks like I'm undercalculating the decay heat quite a bit

[shimwell]

The part I'm most unsure about is this line

https://github.com/jbae11/openmc_activator/blob/5e84aa40044be7cbc288de5df0b3c3f1bd425f4e/openmc_activator.py#L66

I think these fluxes values should be passed in as a user specified argument along with ebins, mg_flux and the others.

I guess Fluxes values are being passed into the other inventory codes, can we use the same numbers for the openmc fluxes

[jbae11]

I'm kinda slow so I still don't think I have a good gasp of what the fluxes should be.. I think it should really just correspond to the material.volume, but please let me know if I'm wrong. I also updated the repo with a mwe and added you so you can try and see if you can make it work :) I think we are super close I'm probably missing something silly