allow overlay material

[shimwell]

This PR allows a material to be used as a tally nuclide bin (virtual overlay material)

multiply_density=False lets a tally respond to a nuclide that is not present in the geometry, and #3771 extended that to void so the response is scored everywhere. That combination is very useful in fusion neutronics: it gives you the dose, damage or reaction rate a detector or a component would see at any point in the model, without perturbing the transport problem by actually putting that material into the geometry.

The gap is that the response is currently per nuclide and microscopic. Real questions are almost always about a material, not a nuclide: absorbed dose in a silicon diode for example. To get that today you have to tally every nuclide of the material separately and recombine the bins by hand, weighted by atom density.

This PR lets you pass the material itself.

What this adds

silicon = openmc.Material()
silicon.add_element('Si', 1.0)
silicon.set_density('g/cm3', 2.33)

tally = openmc.Tally()
tally.filters = [openmc.MeshFilter(mesh)]
tally.scores = ['heating']
tally.nuclides = [silicon]
tally.multiply_density = False

This produces a single bin holding the material's macroscopic response, sum_i N_i * sigma_i, using the overlay material's own atom densities and evaluated everywhere in the geometry including void and regions filled with
something else. For a heating score the result is eV/cm^3 per source particle, so the absorbed dose in Gy is the result multiplied by 1.602e-19 J/eV and divided by the material's mass density in kg/cm^3.

The overlay material does not need to be assigned to any cell. Model picks up materials referenced by tallies and writes them out, so nothing extra is required from the user.

Mixed bins work, for example nuclides = ['Fe56', silicon].

This is a bit of a hack as Tally.nuclides` is now doing two things: it is a list of nuclides, and it can also carry a material. Which is perhaps something we should avoid but the alternative was to add a separate keyword which I think would also not be popular.

Fixes # (issue)

Checklist

• I have performed a self-review of my own code
• I have run clang-format (version 18) on any C++ source files (if applicable)
• I have followed the style guidelines for Python source files (if applicable)
• I have made corresponding changes to the documentation (if applicable)
• I have added tests that prove my fix is effective or that my feature works (if applicable)

[shimwell]

Marking as draft for now to give the F4E and UKAEA people interested in this feature a chance to comment. @alexvalentine94 @Radiation-Transport

[GuySten]

I think that another keyword like materials or virtual_materials is preferable.
We can also restrict to use the transport materials' nuclides or virtual materials but not both of them (to make the implementation a bit easier).
I also think that for virtual materials we can support both cases of multiply_density.

There are also some limitations to virtual materials. I do not think that supporting probability tables is straightforward. We won't be able to support analog estimator like you wrote. But I think we will have a problem with collision estimator especially in void cells.
In void cells a collision estimator won't trigger. And even in non void cells for a tally like heating when the result depend on gamma transport of secondary particles, the result is sensitive to the real material and not the void one so we might want to guard against those cases.

[shimwell]

Thanks @GuySten I shall update and add it via a new kwarg virtual_materials

[paulromano]

Very timely -- I was just working with @j-jeon on a workflow for computing dose in Si so I perfectly understand the need for this! However, I think this would be a lot simpler as a post-processing step on a tally from Python alone. For example, a method that you pass a material and it will multiply the atom densities for each matching nuclide and sum them up in the nuclide dimension. That could be done with a few lines of code rather than the 100s needed here to support it on the C++ side.

[shimwell]

Thanks both for the suggestions.

I should have mentioned this is what we do currently as a python only post process and we are looking for something a little more convenient for the user. Perhaps the additional C++ is too high a cost however, I did get a bit carried away there and like GuySten mention it does not work in various cases. Anyway putting this here for more context

silicon = openmc.Material()
silicon.add_element("Si", 1.0)
silicon.set_density("g/cm3", 2.33)

densities = silicon.get_nuclide_atom_densities()  # {'Si28': 0.04608, ...} atom/b-cm

tally = openmc.Tally(name="si_heating")
tally.filters = [openmc.MeshFilter(mesh)]
tally.scores = ["heating"]
tally.nuclides = list(densities.keys())
tally.multiply_density = False

# run model etc

with openmc.StatePoint(sp_file) as sp:
    t = sp.tallies[tally.id]
    w = np.array([densities[n] for n in t.nuclides])
    macro = (t.mean * w[None, :, None]).sum(axis=1)  # eV/cm3 per source particle
    # eV/cm3 -> Gy, dividing by the overlay material's mass density in kg/cm3
    dose = macro * 1.602176634e-19 / (silicon.get_mass_density() * 1e-3)

[alexvalentine94]

Thanks Jon, this is a very useful feature if implemented that we commonly use in fusion problems. I would agree it is favourable to have this on the input side to aid the reproducibility and lessen the work post processing.