Adding material depletion function

[shimwell]

Description

This PR adds a function that can be used to deplete a bunch of materials with different neutron fluxes, timesteps

Minimal example

mat_ni58 = openmc.Material()
mat_ni58.add_nuclide("Ni58", 0.95)
mat_ni58.set_density("g/cm3", 7.87)
mat_ni58.depletable = True
mat_ni58.volume = 1.0
mat_ni58.temperature = 293.6

depleted_materials = deplete_materials(
    activation_data=[
        {
            "material": mat_ni58,
            "multigroup_flux": np.array([0.5e11] * 42),
            "energy_group_structure": GROUP_STRUCTURES["VITAMIN-J-42"],
            "source_rate": [1e24, 0, 0],
            "timesteps": [10, 10, 10],
        } # could add more entries here
    ],
    timestep_units="s",
    chain_file="chain.xml",
)

@jbae11 and myself have been working on this for a while and I think we have finally brought everything together in the last few weeks.

We have been using the IAEA CoNDERC benchmarks which contain FNS experimental data (decay heat) and FISPACT simulation inputs and outputs for the same irradiation. This benchmark suit is part of the reason that the function API is so flexible as it has irradiation of different materials, with different irradiation schedules with different source rates and different timesteps.

This needed a few additions to PyPact to get the FISPACT inputs parsed nicely. @jbae11 and myself have been working on this since Dec 2023 when we started trying to use the deplete functions https://github.com/fusion-energy/neutronics-workshop/issues/257 It was recently that @jbae11 was able to iron out the main differences between fispact and openmc by making a customized chain file which can be found on the openmc-activator repo I've worked through the openmc-activator converting functions to use pypact and recently adding the same sort of function as I'm adding in this PR. @jbae11 the main differences between this one and the one I added to your repo is that this is more generic as it returns depleted materials and it does a lot of input checking. The one I made for the openmc-activator repo was super useful for testing and produced these comparison plots but I think this one is better in openmc like we discussed. If it is merged in we could make use of it on the openmc-activator repo to continue developing that repo into the main openmc deplete verification and validation repo :rocket: .

Checklist

• [x] I have performed a self-review of my own code
• [x] I have followed the style guidelines for Python source files (if applicable)
• [x] I have made corresponding changes to the documentation (if applicable)
• [x] I have added tests that prove my fix is effective or that my feature works (if applicable)

[lewisgross1296]

Correct me if I'm wrong, but isn't this what IndependentOperator does?

[shimwell]

Correct me if I'm wrong, but isn't this what IndependentOperator does?

It does wrap / build upon that functionality. It uses functionality of MicroXS.from_multigroup_flux, deplete.IndependentOperator and PredictorIntegrator. Just using those functions in python file is however inefficient as that repeatedly loads up the cross sections for each change in material, flux or timesteps. So this PR brings them all together into an openmc.lib call so that it can irradiate multiple materials while only loading the cross sections once which is now fast like the other inventory codes out there.

[shimwell]

This function has been used to produce these comparison results between CoNDERC experimental data and FISPACT simulations https://shimwell.github.io/openmc_activator

[jbae11]

Looks great! Thanks @shimwell. Would it be worth having the option for the user to store the intermediate data like the collapsed cross section?

[shimwell]

taking this PR back to the drawing board, looks like there are some improvements to be made.

I agree @jbae11 we should just return a list of micro_xs instances

I also noticed that I'm currently doing the collapse for all nuclides and not just the ones in the material

I shall tidy up, retest on the @jbae11 openmc-activator repo and modify this PR soon

[shimwell]

To mimic openmc_activator / fispact and other inventory codes the usage of this function is now like this.

The new get_microxs_from_multigroup doesn't add to much compared to the class method MicroXS.from_multigroup_flux. It's main advantages is that it allows one to deplete multiple materials while loading the cross sections once for each nuclide and also returns a list of MicroXS instead of just one. I guess the as MicroXS.from_multigroup_flux is a class method it can't remain a class method and return multiple MicroXS instances. This PR could be seen as a batch version of MicroXS.from_multigroup_flux that works with a multigroup flux instead of performing transport.

Feedback welcome

import openmc
from pathlib import Path
import numpy as np
from openmc.deplete import get_microxs_from_multigroup
from openmc.mgxs import GROUP_STRUCTURES

mat_ni58 = openmc.Material()
mat_ni58.add_nuclide("Ni58", 0.95)
mat_ni58.set_density("g/cm3", 7.87)
mat_ni58.depletable = True
mat_ni58.temperature = 293.6

mat_ni60 = openmc.Material()
mat_ni60.add_nuclide("Ni60", 1.0)
mat_ni60.set_density("g/cm3", 11.34)
mat_ni60.depletable = True
mat_ni60.temperature = 293.6

chain_file_path=Path(__file__).parents[1] / "chain_ni.xml"
multigroup_flux = np.array([0.5e11] * 1102)  # gets normalized in function

all_micro_xs = get_microxs_from_multigroup(  # <----- this PR adds this function here
    materials=openmc.Materials([mat_ni58, mat_ni60]),
    multigroup_fluxes=[multigroup_flux, multigroup_flux],
    energy_group_structures=[
        GROUP_STRUCTURES["UKAEA-1102"],
        "UKAEA-1102",
    ],
    chain_file=chain_file_path,
    reactions=None,
)

for micro_xs, material in zip(all_micro_xs, [mat_ni58, mat_ni60]):
    print('next material')
    material.volume = 1
    operator = openmc.deplete.IndependentOperator(
        materials=openmc.Materials([material]),
        fluxes=[material.volume],
        micros=[micro_xs],
        normalization_mode="source-rate",
        reduce_chain_level=5,
        chain_file=chain_file_path,
    )

    integrator = openmc.deplete.PredictorIntegrator(
        operator=operator,
        timesteps=[1., 10., 10.],
        source_rates=[1e19, 0 ,0],
        timestep_units='s',
    )

    integrator.integrate(path="depletion_results.h5")

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

    for result in results:
        mat = result.get_material(str(material.id))
        print(mat)

updated to use UKAEA-1102 group structure.

[shimwell]

This is an interesting capability. This seems like you are adding utility functions to make O-D transmutation simpler, ala COUPLE-ORIGEN. Is that right?

Yes, just like ALARA, FISPACT, ORIGEN, ACAB etc

I dug into how the multi-group XS collapse is happening, and I see you are using reaction::collapse_rate. From reading the code it doesn't seem to assume any a-priori flux shape. I'm not sure this would be appropriate to use then for relatively sparse energy-group structures (< 1000 groups) like it seems like you are trying to do. Have you looked into if this appropriate collapse treatment?

Yes I think you are right, more groups is better. I have updated the example to use UKAEA-1102 group structure. I have kept the test to a low number to help keep the test quick. We have been able to produce results that are close to fispact and experimental results using 709 groups so I guess the collapse is reasonable.

[yrrepy]

I think this kind of functionality is great. Allowing for simple, straightforward activation. Few people know that OpenMC has activation in it

[shimwell]

I have updated this PR to take advances of the recent changes to the chain file loading.

I have also profiled the example simulation when running the above example passing in a chain as a path/string an passing it in as a preloaded Chain object

Time as string = 18 seconds Time as Chain = 9 seconds

Full details, images, reproducible scripts here https://github.com/shimwell/openmc_profile

@jbae11 and @eepeterson might be interested in seeing this speed improvement.

[eepeterson]

I like this work a lot. I do have a few comments on the interface though that I'd like to get to over the next few days sometime. The primary one is that I think a functionality like this should live on the Material and Materials classes as a deplete method that takes in the flux, energy group structure and operator and integrator kwargs. Some of these arguments are redundant and/or unnecessary in a simplified activation calculation like this and could be greatly cleaned up for the user which might be nice to see included in this PR. I could also see a benefit of not writing out a depletion_results.h5 file and just passing back (in the case this is called on a single material) a list of materials the same length as the time step vector or +1 for the original composition whichever you like. For now, (barring additional refactoring of the depletion module) you could call the integrate function from a temp directory and just return the depleted materials list. Also, as a todo we should note that in the future we will probably want to add options for different in group flux treatments which I think FISPACT has and gets to @MicahGale's point.

All of the above is not intended to be instead of the current get_microxs_from_multigroup method, but on top of it for a cleaner workflow.

[shimwell]

@eepeterson I've added the material.deplete as suggested, this is something I had planned to do in a follow up PR but happy to add it here.

the minimal usage is

    pristine_material = openmc.Material()
    pristine_material.add_nuclide("Ni58", 0.95)
    pristine_material.set_density("g/cm3", 7.87)
    pristine_material.depletable = True  # could automate this in the function
    pristine_material.temperature = 293.6
    pristine_material.volume = 1.

    mg_flux = [0.5e11] * 42

    depleted_material = pristine_material.deplete(
        multigroup_flux=mg_flux,
        energy_group_structures="VITAMIN-J-42",
        timesteps=[10,10],
        source_rates=[1e19, 0.0],
    )

i shall add openmc.Materials.deplete soon

[shimwell]

Just added the materials.deplete. Could do with some good tests and I need to check it through for typos but I think we are closing in on a nice user API

I think I can reduce the code a bit with another round. Will do this soon

[shimwell]

I've tidied up the PR and added a few more tests to check the depleted materials are as expected.

Minimal usage for depleting and openmc.Material or openmc.Materials is

import openmc

openmc.config['chain_file'] = '/home/jon/nuclear_data/chain-endf-b8.0.xml'
openmc.config['cross_sections'] = '/home/jon/nuclear_data/cross_sections.xml'

mg_flux = [0.5e11] * 42

pristine_material_1 = openmc.Material(material_id=1)
pristine_material_1.add_nuclide("Ni58", 1.)
pristine_material_1.set_density("g/cm3", 7.87)
pristine_material_1.temperature = 293.6
pristine_material_1.volume = 1.

# deplete a single material
depleted_material = pristine_material_1.deplete(
    multigroup_flux=mg_flux,
    energy_group_structure="VITAMIN-J-42",
    timesteps=[100,100],
    source_rates=[1e19, 0.0],
    timestep_units='d',
)

pristine_material_2 = openmc.Material(material_id=2)
pristine_material_2.add_nuclide("Ni60", 1.)
pristine_material_2.set_density("g/cm3", 7.87)
pristine_material_2.temperature = 293.6
pristine_material_2.volume = 1.

# deplete multiple materials material
pristine_materials = openmc.Materials([pristine_material_1, pristine_material_2])

depleted_materials = pristine_materials.deplete(
    multigroup_fluxes=[mg_flux, mg_flux],
    energy_group_structures=["VITAMIN-J-42", "VITAMIN-J-42"],
    timesteps=[10, 76000],
    source_rates=[1e19, 0.0],
    timestep_units='a',
)

Keen to hear what people think.

[shimwell]

@shimwell Thanks for this PR. I really like the concept but I do see that the get_microxs_from_multigroup has quite a bit of duplicated code from the from_multigroup_flux classmethod. I have an idea for how to eliminate this duplication with a new context manager, but it might require expanding the scope of this PR, so I think it's probably better left as a separate PR. Are you OK with me submitting a new PR separately and then we can come back to this PR and refactor?

Yes please I would be keen to make use of a new context manager. I shall look out for your PR and refactor this branch after your context manager gets merged in. Many thanks

[shimwell]

I have now made use of the temporary session context manager. This also allowed me to remove a lot of code introduced earlier in the PR.

[shimwell]

Thanks Paul that is some really nice tidying up. I particularly appreciate the speed increase and the nuclide part. Yes very happy to merge. (Auto Merge now enabled)