Python API not reading nuclear data ```interpolation``` and ```breakpoint``` correctly

[shimwell]

Bug Description

Reading in the h5 nuclear data with openmc appears to result in default interpolation and breakpoint values such as [2]

Steps to Reproduce

change the dir in this script and run the script it should print out different interpolation and breakpoint values when the file is read with the openmc python api and h5py directly.

import openmc.data
import h5py

# Load a sample HDF5 file and product
filename = "/home/jon/nuclear_data/tendl-2021-hdf5/tendl-2021-hdf5/Cr52.h5"
lib = openmc.data.IncidentNeutron.from_hdf5(filename)

# Pick a product with tabulated distributions
for rx in lib.reactions.values():
    for prod in rx.products:
        # Look for Tabulated1D in applicability or yield
        for app in getattr(prod, "applicability", []):
            if hasattr(app, "breakpoints") and hasattr(app, "interpolation"):
                print("App breakpoints:", app.breakpoints)
                print("App interpolation:", app.interpolation)
        y = getattr(prod, "yield_", None)
        if hasattr(y, "breakpoints") and hasattr(y, "interpolation"):
            print("Yield breakpoints:", y.breakpoints)
            print("Yield interpolation:", y.interpolation)
        # Look for Tabulated1D in distributions
        for dist in getattr(prod, "distribution", []):
            if hasattr(dist, "breakpoints") and hasattr(dist, "interpolation"):
                print("Dist breakpoints:", dist.breakpoints)
                print("Dist interpolation:", dist.interpolation)

# Now print the same info using h5py
print("\n--- h5py direct read ---")
with h5py.File(filename, "r") as f:
    cr52 = f["Cr52"]
    reactions = cr52["reactions"]
    for rx_name in reactions:
        rx_group = reactions[rx_name]
        for prod_name in rx_group:
            if not prod_name.startswith("product_"):
                continue
            prod_group = rx_group[prod_name]
            # Yield
            if "yield" in prod_group:
                y_ds = prod_group["yield"]
                print(f"{rx_name} {prod_name} yield breakpoints:", y_ds.attrs.get("breakpoints", None))
                print(f"{rx_name} {prod_name} yield interpolation:", y_ds.attrs.get("interpolation", None))
            # Distribution(s)
            for dist_key in prod_group:
                if dist_key.startswith("distribution_"):
                    dist_group = prod_group[dist_key]
                    # Applicability
                    if "applicability" in dist_group:
                        app_ds = dist_group["applicability"]
                        print(f"{rx_name} {prod_name} {dist_key} applicability breakpoints:", app_ds.attrs.get("breakpoints", None))
                        print(f"{rx_name} {prod_name} {dist_key} applicability interpolation:", app_ds.attrs.get("interpolation", None))
                    # Angle/mu distribution
                    if "angle" in dist_group and "mu" in dist_group["angle"]:
                        mu_ds = dist_group["angle"]["mu"]
                        print(f"{rx_name} {prod_name} {dist_key} angle/mu interpolation:", mu_ds.attrs.get("interpolation", None))

Environment

ubuntu , python 3.13, openmc develop branch

[GuySten]

This script print a large amount of numbers. Could you provide a specific reaction and reaction product where there is a problem. It will help if you could also provide the script output for that reaction and reaction product combination using python api and h5py.

[shimwell]

Thanks GuySten

MT 70 products appear to have different break points, others are also different but this script focuses on mt 70

import openmc.data
import h5py

# Load a sample HDF5 file and product
filename = "/home/jon/nuclear_data/tendl-2021-hdf5/tendl-2021-hdf5/Cr52.h5"
lib = openmc.data.IncidentNeutron.from_hdf5(filename)

reaction_mt = 70

# Pick a product with tabulated distributions
for rx in lib.reactions.values():
    for prod in rx.products:
        if rx.mt == reaction_mt:
            # Look for Tabulated1D in applicability or yield
            for app in getattr(prod, "applicability", []):
                if hasattr(app, "breakpoints") and hasattr(app, "interpolation"):
                    print("App breakpoints:", app.breakpoints)
                    print("App interpolation:", app.interpolation)
            y = getattr(prod, "yield_", None)
            if hasattr(y, "breakpoints") and hasattr(y, "interpolation"):
                print("Yield breakpoints:", y.breakpoints)
                print("Yield interpolation:", y.interpolation)
            # Look for Tabulated1D in distributions
            for dist in getattr(prod, "distribution", []):
                if hasattr(dist, "breakpoints") and hasattr(dist, "interpolation"):
                    print("Dist breakpoints:", dist.breakpoints)
                    print("Dist interpolation:", dist.interpolation)

# Now print the same info using h5py
print("\n--- h5py direct read ---")
with h5py.File(filename, "r") as f:
    cr52 = f["Cr52"]
    reactions = cr52["reactions"]
    for rx_name in reactions:
        if rx_name == f"reaction_0{reaction_mt}":
            rx_group = reactions[rx_name]
            for prod_name in rx_group:
                if not prod_name.startswith("product_"):
                    continue
                prod_group = rx_group[prod_name]
                # Yield
                if "yield" in prod_group:
                    y_ds = prod_group["yield"]
                    # print(f"{rx_name} {prod_name} yield breakpoints:", y_ds.attrs.get("breakpoints", None))
                    print(f"{rx_name} {prod_name} yield interpolation:", y_ds.attrs.get("interpolation", None))
                # Distribution(s)
                for dist_key in prod_group:
                    if dist_key.startswith("distribution_"):
                        dist_group = prod_group[dist_key]
                        # Applicability
                        if "applicability" in dist_group:
                            app_ds = dist_group["applicability"]
                            print(f"{rx_name} {prod_name} {dist_key} applicability breakpoints:", app_ds.attrs.get("breakpoints", None))
                            print(f"{rx_name} {prod_name} {dist_key} applicability interpolation:", app_ds.attrs.get("interpolation", None))
                        # Angle/mu distribution
                        if "angle" in dist_group and "mu" in dist_group["angle"]:
                            mu_ds = dist_group["angle"]["mu"]
                            print(f"{rx_name} {prod_name} {dist_key} angle/mu interpolation:", mu_ds.attrs.get("interpolation", None))

for Cr52.h5 Tendl 2021 outputs this

App breakpoints: [2]
App interpolation: [2]
Yield breakpoints: [2]
Yield interpolation: [2]
Yield breakpoints: [2]
Yield interpolation: [2]
Yield breakpoints: [2]
Yield interpolation: [2]

--- h5py direct read ---
reaction_070 product_0 yield interpolation: None
reaction_070 product_0 distribution_0 applicability breakpoints: [2]
reaction_070 product_0 distribution_0 applicability interpolation: [2]
reaction_070 product_0 distribution_0 angle/mu interpolation: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]
reaction_070 product_1 yield interpolation: [2]
reaction_070 product_1 distribution_0 angle/mu interpolation: [1 1]
reaction_070 product_2 yield interpolation: [2]
reaction_070 product_2 distribution_0 angle/mu interpolation: [1 1]
reaction_070 product_3 yield interpolation: [2]
reaction_070 product_3 distribution_0 angle/mu interpolation: [1 1]

the python api shows lin lin [2] in all cases where the hdf5 data contains histogram [1, 1] in some cases

[GuySten]

I think the problem is in your script: The script does not know how to parse angle interpolation.

You can see that in the output there are no lines with either Dist breakpoints or Dist interpolation.

[shimwell]

Sorry my script is a mess, I think I've found something in the openmc source code

the Reation.from_hdf5 function calls Tabulated1D with no arguments for breakpoints or interpolation

https://github.com/openmc-dev/openmc/blob/f70febb05f5fc32925d08707db48ff2c33dfe493/openmc/data/reaction.py#L975C21-L975C78

So both breakpoints or interpolation default to None and then we get default of [2] linear linear interpolation

https://github.com/openmc-dev/openmc/blob/f70febb05f5fc32925d08707db48ff2c33dfe493/openmc/data/function.py#L147-L148

I think this is not correct to allow defaults here

[GuySten]

I think it is okay because the cross sections that are stored in openmc hdf5 libraries always have linear-linear interpolation.

EDIT: One of the responsibilities of NJOY2016 is to convert cross sections to linear-linear representation.

[shimwell]

I agree on the cross sections. However this is for the reaction products that have their own distributions.

I think the reaction products are missing their interpolation values when being read from hdf5 in the python api. When writing products in the python api that looks fine and when reading products in the c++ code that also looks fine.

[GuySten]

I think the reaction products are missing their interpolation values when being read from hdf5 in the python api. When writing products in the python api that looks fine and when reading products in the c++ code that also looks fine.

I did not find evidence to support that. I could look into it if you find some evidence.

[shimwell]

Thanks for your help with this, much appreachiated

I include a screenshot of h5dump of Cr52.h5 where we can see the interpolation is [1,1] for reaction 70.

When reading in this nuclear data with the openmc python api it gives [2]

This is the same reaction as my earlier script was printing out interpolation values for.

[GuySten]

But how do you know that the python api does not read this correctly?