Plotting flux spectrum after depletion in openmc

Here is how I can view the neutron flux at the output. There is an application in Jupyter notebook that is made like this. @paulromano @Shimwell @pshriwise

sp_path = model.run(output=False)

How can I change this for the following integrator?

# Instantiate a Settings object, set all runtime parameters, and export to XML
settings = openmc.Settings()
settings.batches = 500
settings.inactive = 20
settings.particles = 500

# Create an initial uniform spatial source distribution over fissionable zones
settings.source = openmc.source.Source(space=openmc.stats.Point())

tally = dict()
tally['heating'] = openmc.Tally(name="heating")
tally['heating'].scores.append('heating-local')
e_min, e_max = 1e-5, 20e6
groups = 500
energies = np.logspace(log10(e_min), log10(e_max), groups + 1)
energy_filter = openmc.EnergyFilter(energies)
particle_filter = openmc.ParticleFilter(['neutron'])
cell_filter = openmc.MaterialFilter([uo2])
tally['flux']= openmc.Tally(name="flux")
tally['flux'].filters = [energy_filter, particle_filter] 
tally['flux'].scores = ['flux']
mesh = openmc.RegularMesh()
mesh.dimension = [500, 500, 1]
mesh.lower_left = [-100, -100, 50]
mesh.upper_right = [100, 100, 200]
mesh_filter = openmc.MeshFilter(mesh)
tally['mesh'] = openmc.Tally(name="Mesh")
tally['mesh'].scores = ['flux','absorption','fission']
tally['mesh'].filters = [mesh_filter]
tally['mesh'].filters.append(particle_filter)
tallies = openmc.Tallies(tally.values())

###############################################################################
#                   Initialize and run depletion calculation
###############################################################################

with open('/depletion/serpent_fissq.json', 'r') as f:
    serpent_fission_q = json.load(f)

model = openmc.Model(geometry=geom, settings=settings, tallies=tallies)
model.export_to_xml()
chain_file = '/depletion-comparison/data/depletion/chain.xml'
op = openmc.deplete.Operator(model, chain_file,
    fission_q=serpent_fission_q,
    fission_yield_mode="average")

# cumulative steps in MWd/kg
burnup_cum = np.array([0.1, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.5, 15.0, 17.5, 20.0, 22.5, 25.0, 27.5, 30.0, 32.5, 35.0, 37.5, 40.0, 42.5, 45.0])
burnup = np.diff(burnup_cum, prepend=0.0)
power = 131040  

integrator = openmc.deplete.PredictorIntegrator(op, burnup, power, timestep_units='MWd/kg')
integrator.integrate()

hi prof electron, the specified tally will be reported on the statepoint file of each depletion steps, so to read your tally output like from the second depletion steps, you could read your statepoint_n2.h5, and so on. so as long as your tally has been added to the model for your transport operator being used at the integrator, then it will be good.
model = openmc.Model(geometry, materials, settings, tallies)
op = openmc.deplete.Operator(model, …
you can try it with small particle history just to make sure that everything is as you want.

2 Likes

Hi, sorry to disturb you. I tryed add the tally into my model, but the ‘statepoint_nx.h5’ didn’t exest after whole 500 cycles calculation. There are still ‘openmc_simulation_nx.h5’ and ‘statepoin.500.h5’ only.

my code was:
model = openmc.Model(geometry=geom, settings=settings, tallies=tallies_file)
operator = openmc.deplete.CoupledOperator(model, ‘/home/ads/vir-jiang/openmc/chain_endfb80_sfr.xml’)

maybe I should add the materials=xx to the model? :face_with_monocle:

Hi Jiang
Sorry for the confusion, yes, the openmc statepoint file for each depletion (burnup) step was named openmc_simulation_nx.h5. So it’s not statepoint_nx.h5. I just want to said that the statepoint file was produced for each step with its unique name so you can use it later.
About adding materials to the openmc.model, I think it wasn’t necessary because openmc model will regenerate the material.xml file later.
Hope it helps you

1 Like

Thanks a lot. :wink:
Yes, it’s not the reason of “material”, I type "sp1 = openmc.StatePoint(‘openmc_simulation_n1.h5’) " and it worked.

1 Like