Possible Issue in Decay-Only Simulation with Fixed-Source Transmutation

User Support
1

Hi OpenMC experts,

I’ve encountered a problem when using OpenMC to study a Fixed-Source Transmutation problem and the associated decay photon spectrum.

When trying to simulate a decay-only scenario during the whole activation process by setting the source rate to zero, negative nuclide densities appear. I also noticed that this issue occurs particularly when only the nuclide H1 (e.g., in water) is defined.

Other materials seem to work fine so far.

WARNING: nuclide Cd107 in material2 is negative (density = -1.9449206001183648e-19 atom/b-cm)
WARNING: nuclide Sm146 in material2 is negative (density = -1.6800623999344763e-21 atom/b-cm)

Moreover, this leads to an error when using the get_decay_photon_energy() function. In a decay-only water problem(2 time steps, 1 day/step), we shouldn’t get any photon intensity, but the function still returns incorrect results.

Just wanted to share this situation — any suggestions or insights would be greatly appreciated.

Thanks!

import openmc
import openmc.deplete
import math

water = openmc.Material(2, name="water")
water.depletable = True
water.volume = 20 * 20 * 20  # cm³
water.add_nuclide('H1', 2.0)
water.add_nuclide('O16', 1.0)
water.set_density('g/cm3', 1.0)

materials = openmc.Materials([water])

x_min = openmc.XPlane(x0=-10,  boundary_type='vacuum')
x_max = openmc.XPlane(x0=10,  boundary_type='vacuum')
y_min = openmc.YPlane(y0=-10,  boundary_type='vacuum')
y_max = openmc.YPlane(y0=10,  boundary_type='vacuum')
z_min = openmc.ZPlane(z0=-10,  boundary_type='vacuum')
z_max = openmc.ZPlane(z0=10,  boundary_type='vacuum')

cube = +x_min & -x_max & +y_min & -y_max & +z_min & -z_max
cube = openmc.Cell(region=cube, fill=water)

universe = openmc.Universe(cells=[cube])
geometry = openmc.Geometry(universe)

source = openmc.IndependentSource()
source.space = openmc.stats.Point((0, 0, 0))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Discrete([14e6], [1.0])
source.particle = 'neutron'

settings = openmc.Settings()
settings.batches = 2
settings.inactive = 0
settings.particles = 100000
settings.source = source
settings.run_mode = 'fixed source'

model = openmc.Model(geometry, materials, settings)
chain_file = "./chain_endfb80_pwr.xml"
operator = openmc.deplete.CoupledOperator(model, chain_file, normalization_mode='source-rate')

test_timesteps = [1 * 24 *3600, 1 * 24 *3600]
source_rates = [0,0]
integrator = openmc.deplete.PredictorIntegrator(operator, test_timesteps, source_rates=source_rates)
integrator.integrate()
2

The CRAM solver does this. Those nuclides are filtered out before ending up in the solution. See code block below.

Post more about your get_decay_photon_energy error if you want assistance with that.

openmc/openmc/deplete/independent_operator.py

    def _update_materials(self):
        """Updates material compositions in OpenMC on all processes."""

        for rank in range(comm.size):
            number_i = comm.bcast(self.number, root=rank)

            for mat in number_i.materials:
                nuclides = []
                densities = []
                for nuc in number_i.nuclides:
                    if nuc in self.nuclides_with_data:
                        val = 1.0e-24 * number_i.get_atom_density(mat, nuc)

                        # If nuclide is zero, do not add to the problem.
                        if val > 0.0:
                            if self.round_number:
                                val_magnitude = np.floor(np.log10(val))
                                val_scaled = val / 10**val_magnitude
                                val_round = round(val_scaled, 8)

                                val = val_round * 10**val_magnitude

                            nuclides.append(nuc)
                            densities.append(val)
                        else:
                            # Only output warnings if values are significantly
                            # negative. CRAM does not guarantee positive
                            # values.
                            if val < -1.0e-21:
                                print(f'WARNING: nuclide {nuc} in material'
                                      f'{mat} is negative (density = {val}'

                                      ' atom/b-cm)')
                            number_i[mat, nuc] = 0.0
3

Thanks for your kind and prompt reply. :blush:

My get_decay_photon_energy() function is written as follows:

openmc.config['chain_file'] = "chain_endfb80_pwr.xml"

results = openmc.deplete.Results("depletion_results.h5")
step = results[-1]

activated_mat = step.get_material('2')

photon_spectrum = activated_mat.get_decay_photon_energy(clip_tolerance=1e-6,units="Bq")

energies_eV = photon_spectrum.x
intensities = photon_spectrum.p  

#eV to MeV
energies = np.array(energies_eV) * 1e-6

for e, I in zip(energies, intensities):
    print(f"Energy: {e:.6f} MeV, Intensity: {I:.6e} Bq")

print(f"Total intensity: {np.sum(intensities):.6e}")
print("Material id:", activated_mat.id)

In the case I mentioned earlier, the model is a cube filled with water, and I performed two steps of decay-only calculation.
Since there shouldn’t be any decay reactions in this case, the get_decay_photon_energy() function should ideally return no output.
However, it still produces the following results, which are clearly incorrect.

Screenshot 2025-10-29 at 02.18.19
Screenshot 2025-10-29 at 02.18.19802×362 10.9 KB

I’ve checked the user guide, but I’m not sure if there’s a way to filter out negative density values to avoid such wrong outputs from get_decay_photon_energy().

I discovered this issue while studying activation scenarios such as:

sources_rate = [1e10, 0, 1e10, 0]

This setup is quite common in activation or SDR calculations.
For cells filled with steel, the results are consistent, but for water, the issue appears since the 0 sources rate as described.

4

You might try using independent operator.

The Bq are very low.
Another particularity is that depletion populates other nuclides that are not part of the problem.

Probably it created these spurious nuclides. Inspect the depletion_results.h5 for these nuclides that don’t belong.

They need to be filtered out from the results.