Discrepant results using surf_source_write / surf_source_read

User Support
1

Hi,

I am doing a very simple simulation with a neutron isotropic point source at the centre of nested water spheres to compare the current and flux directly tallied with those obtained, first writing to a surface source file (surf_source_write), then reading from this file (surf_source_read). Shouldn’t I get the same result in both cases? Maybe I am doing something stupid but I don’t understand why I obtain completely different results.

The file for the simulation is really simple:

import openmc

iwrite = True # write surface source
#iwrite = False # read surface source

water = openmc.Material()
water.add_element(‘H’,2,‘ao’)
water.add_element(‘O’,1,‘ao’)
water.set_density(‘g/cm3’,1.)
materials_file = openmc.Materials([water])
materials_file.export_to_xml()

s1 = openmc.Sphere(surface_id=1,r= 5.)
s2 = openmc.Sphere(surface_id=2,r=10.)
s3 = openmc.Sphere(surface_id=3,r=15.)
s4 = openmc.Sphere(surface_id=4,r=20.)
s5 = openmc.Sphere(surface_id=5,r=25.,boundary_type=‘vacuum’)
cell1 = openmc.Cell(fill=water, region= -s1)
cell2 = openmc.Cell(fill=water, region= +s1 & -s2)
cell3 = openmc.Cell(fill=water, region= +s2 & -s3)
cell4 = openmc.Cell(fill=water, region= +s3 & -s4)
cell5 = openmc.Cell(fill=water, region= +s4 & -s5)
geometry_file = openmc.Geometry([cell1,cell2,cell3,cell4,cell5])
geometry_file.export_to_xml()

settings_file = openmc.Settings()
settings_file.batches = 30
settings_file.particles = 500
settings_file.run_mode = ‘fixed source’
if (iwrite):
source = openmc.IndependentSource()
source.space = openmc.stats.Point(xyz=(0,0,0))
source.angle = openmc.stats.Isotropic()
settings_file.source = source
settings_file.surf_source_write = {‘surface_ids’: [2],‘max_particles’: 15000}
else:
settings_file.surf_source_read = {‘path’: ‘surface_source.h5’}
settings_file.export_to_xml()

surf_filter = openmc.SurfaceFilter([s1,s2,s3,s4,s5])
cell_filter = openmc.CellFilter([cell1,cell2,cell3,cell4,cell5])
tallies_file = openmc.Tallies()
tally = openmc.Tally()
tally.filters = [surf_filter]
tally.scores = [‘current’]
tallies_file.append(tally)
tally = openmc.Tally()
tally.filters = [cell_filter]
tally.scores = [‘flux’]
tallies_file.append(tally)
tallies_file.export_to_xml()

openmc.run()

2

Hi Josem,

in the SSW, the geometry outside of the surface source surface (2) should be voided.

I would compare a full, non SSW/SSR run with the SSR run.
Only compare tallies for cells/surfaces beyond the writing/reading surface.

The full run tallies will be normalized to it’s particles.
The SSR run tallies will need to be renormalized to that, so:
tally_full=tally_ssr * (SSR_particles) / (SSW_particles)

1 Like
3

Hi,
Thank you very much! Your explanation seems now pretty obvious but I’ve been thinking about the question all the morning in vain!!

4

Using HDFView I can see the number of particles stored in the surface_source.h5 file. Is there any command to directly obtain this number ?

5

I mostly use the new MCPL format files and then use MCPL itself to inspect my phase-space.

I usually run my SSR some multiple or some fraction of what is on the phase-space.
But if I recall correctly, you don’t actually need to know how many are on the phase-space it for your normalization.
But I might not be recalling well, this stuff can get confusing.

You can use the following Python code to inspect your surf_source.h5

import h5py
import os
import numpy as np
import openmc

source_file = h5py.File('surface_source.h5')
particles = source_file['source_bank'][()]

print(len(particles))      # number of particles on file
1 Like