I am currently running some tests in order to check the correct transfer of mass from one area to another.
I have a sphere of 239Pu used solely for transport, that won’t be of interest here.
The main point is a cloud made of 4 gases out of the geometry, followed by 4 reservoirs.
The 4 gases are extracted from the cloud towards the 1st reservoir, then 3 are extracted from reservoir 1 towards the next etc.
When I extract nucleide masses and plot them, I observe that some gas masses are reset to their original value in the cloud, with no consequence on the correct distribution in the reservoirs.
Has anyone else seen such behaviour ?
Input and post-processing python files are provided below.
I asked one collegue to do the same with his different config and he observed such things as well. I also tried with other fuels and other gases and the same behaviour occurs, always with the first material, whether it is included in the geometry or not.
My thanks in advance to any one contributing to this post 
The input :
import openmc
import pandas as pd
import numpy as np
import openmc.deplete
chain_file = "" #a chain file
xs_file = "" #a xs file, I used data from ENDF-B7.1 for this test
openmc.config["cross_sections"] = xs_file
openmc.config["chain_file"] = chain_file
temperature = 900
# Geometry
radius_core = 4.95
core = openmc.Sphere(r=radius_core, boundary_type="vacuum")
# Cells
incore = openmc.Cell(name="incore", region=-core)
outcore = openmc.Cell(name="outcore", region=+core)
mat_fuel = openmc.Material(name="fuel", temperature=temperature)
mat_fuel.add_components({
"Pu239": 1.0,
})
mat_fuel.set_density("g/cm3", 19.8)
mat_fuel.volume = (4/3)*np.pi*(radius_core**3)
mat_fuel.depletable = True
incore.fill = mat_fuel
mat_cloud = openmc.Material(name="cloud", temperature=temperature)
mat_cloud.add_components({
"He4": 0.25,
"Ar40": 0.25,
"Kr84": 0.25,
"Xe132": 0.25
})
radius_cloud = 10
mat_cloud.set_density("g/cm3", 1e-3)
mat_cloud.volume = (4/3)*np.pi*(radius_cloud**3)
mat_cloud.depletable = True
# Universe
root_universe = openmc.Universe(cells=(incore, outcore))
geometry = openmc.Geometry()
geometry.root_universe = root_universe
# Settings
# Create a point source
point = openmc.stats.Point((0, 0, 0))
source = openmc.Source(space=point)
settings = openmc.Settings()
settings.source = source
settings.batches = 400
settings.inactive = 100
settings.particles = 5000
settings.output = {"path":"."}
# Depletion
Reservoir1 = openmc.Material(name='Reservoir1')
Reservoir1.add_nuclide("Fe56", 1) # random stable element to avoid empty mat
Reservoir1.set_density("g/cm3", 1e-32)
Reservoir1.volume = 1e6
Reservoir1.depletable = True
Reservoir2 = openmc.Material(name='Reservoir2')
Reservoir2.add_nuclide("Fe56", 1)
Reservoir2.set_density("g/cm3", 1e-32)
Reservoir2.volume = 1e6
Reservoir2.depletable = True
Reservoir3 = openmc.Material(name='Reservoir3')
Reservoir3.add_nuclide("Fe56", 1)
Reservoir3.set_density("g/cm3", 1e-32)
Reservoir3.volume = 1e6
Reservoir3.depletable = True
Reservoir4 = openmc.Material(name='Reservoir4')
Reservoir4.add_nuclide("Fe56", 1)
Reservoir4.set_density("g/cm3", 1e-32)
Reservoir4.volume = 1e6
Reservoir4.depletable = True
materials = openmc.Materials([mat_fuel, mat_cloud, Reservoir1, Reservoir2, Reservoir3, Reservoir4])
materials.export_to_xml()
geometry.export_to_xml()
settings.export_to_xml()
times = np.array([1/24, 3/24, 6/24, 12/24, 1, 2, 4, 8, 16, 32, 64, 128, 256, 384, 512, 640, 768, 896, 1024])
times2 = times*1
for i in range(1, len(times2)):
times2[i] = times[i] - times[i-1]
power = 10e6
towards_Res1 = ["He", "Ar", "Kr", "Xe"]
towards_Res2 = ["He", "Ar", "Kr"]
towards_Res3 = ["He", "Ar"]
towards_Res4 = ["He"]
model = openmc.model.Model(geometry, materials, settings)
operator = openmc.deplete.CoupledOperator(model, chain_file=chain_file, normalization_mode='fission-q', fission_yield_opts={"energy":500e3})
integ = openmc.deplete.PredictorIntegrator(operator, times2, power, timestep_units='d')
integ.add_transfer_rate(mat_cloud, towards_Res1, 1/(1*3600), destination_material=Reservoir1)
integ.add_transfer_rate(Reservoir1, towards_Res2, 1/(2*3600), destination_material=Reservoir2)
integ.add_transfer_rate(Reservoir2, towards_Res3, 1/(3*3600), destination_material=Reservoir3)
integ.add_transfer_rate(Reservoir3, towards_Res4, 1/(4*3600), destination_material=Reservoir4)
integ.integrate()
The post-processing :
import openmc
import openmc.deplete
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
chain_file = "" # same file as input
xs_file = "" # same file as input
openmc.config["cross_sections"] = xs_file
openmc.config["chain_file"] = chain_file
f_dep = "depletion_results.h5"
results = openmc.deplete.Results(f_dep)
nuclei = ["He4", "Ar40", "Kr84", "Xe132"]
colors = ["red", "orange", "cyan", "blue"]
areas = ["Cloud", "Reservoir 1", "Reservoir 2", "Reservoir 3", "Reservoir 4"]
times = np.array([0, 1/24, 3/24, 6/24, 12/24, 1, 2, 4, 8, 16, 32, 64, 128 , 256, 384, 512, 640, 768, 896, 1024])
dico_results = {}
for area in areas:
dico_results[area] = {}
for nuc in nuclei:
dico_results[area][nuc] = np.zeros((len(times)))
for id_mat in range(1, len(areas)+1): #core is 0
area = areas[id_mat-1]
for i in range(len(times)):
mat = results.export_to_materials(i)[id_mat]
nucl = mat.get_nuclides()
# print(dico_results[area].keys())
# print(nucl)
for n in range(len(nucl)):
try: dico_results[area][nucl[n]][i] = mat.get_mass(nucl[n])
except: print(f"{nucl[n]} : {mat.get_mass(nucl[n])}") # ignore any element that can appear with superlow value, including the very low amount of iron
# print(dico_results.keys())
for area in areas:
for nuc in nuclei:
print(dico_results[area][nuc])
fig = plt.figure(figsize=(4*len(areas), 4))
gs = GridSpec(1, len(areas), fig)
axs=[]
for i in range(len(areas)):
axs.append(fig.add_subplot(gs[i]))
for i in range(len(areas)):
for j in range(len(nuclei)):
vect = dico_results[areas[i]][nuclei[j]]
axs[i].scatter(times, vect, label=nuclei[j], color=colors[j])
axs[i].set_xscale("log")
axs[i].set_xlim([times[1]*0.5, times[-1]*1.2])
axs[i].set_xlabel("Time [d]")
axs[i].set_ylabel("Mass [g]")
axs[i].set_title(areas[i])
axs[0].legend()
fig.tight_layout()
fig.savefig("simple_distribution_in_areas.png", dpi=300)