# Discrepant results using surf\_source\_write / surf\_source\_read

**URL:** <https://openmc.discourse.group/t/discrepant-results-using-surf-source-write-surf-source-read/3412>\
**Category:** User Support\
**Created:** [September 11, 2023, 1:57pm UTC](https://openmc.discourse.group/t/discrepant-results-using-surf-source-write-surf-source-read/3412 "2023-09-11T13:57:05Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![josem](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/josem/32/4096_2.png) [@josem](https://openmc.discourse.group/u/josem)\
**Post date:** [September 11, 2023, 1:57pm UTC](https://openmc.discourse.group/t/discrepant-results-using-surf-source-write-surf-source-read/3412/1 "2023-09-11T13:57:05Z")

</div>

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()

---

<div class="post-metadata">

**Author:** ![yrrepy](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/yrrepy/32/2365_2.png) [@yrrepy](https://openmc.discourse.group/u/yrrepy)\
**Post date:** [September 11, 2023, 4:01pm UTC](https://openmc.discourse.group/t/discrepant-results-using-surf-source-write-surf-source-read/3412/2 "2023-09-11T16:01:21Z")

</div>

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)

---

<div class="post-metadata">

**Author:** ![josem](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/josem/32/4096_2.png) [@josem](https://openmc.discourse.group/u/josem)\
**Post date:** [September 11, 2023, 8:03pm UTC](https://openmc.discourse.group/t/discrepant-results-using-surf-source-write-surf-source-read/3412/3 "2023-09-11T20:03:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![josem](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/josem/32/4096_2.png) [@josem](https://openmc.discourse.group/u/josem)\
**Post date:** [September 11, 2023, 11:01pm UTC](https://openmc.discourse.group/t/discrepant-results-using-surf-source-write-surf-source-read/3412/4 "2023-09-11T23:01:26Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![yrrepy](https://yyz2.discourse-cdn.com/free1/user_avatar/openmc.discourse.group/yrrepy/32/2365_2.png) [@yrrepy](https://openmc.discourse.group/u/yrrepy)\
**Post date:** [September 12, 2023, 12:01am UTC](https://openmc.discourse.group/t/discrepant-results-using-surf-source-write-surf-source-read/3412/5 "2023-09-12T00:01:11Z")

</div>

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

```auto
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

```