Spherical photon source - discrepancies with MCNP and Serpent

Hi!

I’m interested in using OpenMC to determine the dose rate outside a spherical Am241 source with a radius of 1 cm. I set up the spherical spatial distribution of the source as follows:

sphere_radius = 1
r = openmc.stats.PowerLaw(0, sphere_radius, 1)
theta = openmc.stats.Uniform(0, pi)
phi = openmc.stats.Uniform(0, 2*pi)
spatial_dist = openmc.stats.SphericalIndependent(r, theta, phi)
source.space = spatial_dist

Originally, I used a source energy of 60 keV (the main line from Am241), but got results that seemed a bit to low when comparing with references. To cross-check the method, I increased the source energy to 1 MeV (still sampling from the spherical distribution), and implemented the same geometry and source in Serpent and MCNP for benchmarking. I tally the photon current going out of the 1-cm sphere. The results are shown below (the current across the sphere surface per initial source particle):

Spherical source:
OpenMC: 0.5898 (current in 1 MeV bin: 0.4694)
Serpent: 0.6218 (current in 1 MeV bin: 0.5011)
MCNP: 0.6213 (current in 1 MeV bin: 0.5012)

Clearly, OpenMC yields both a lower total current and a lower fraction of the uncollided 1-MeV photons passing through the sphere surface than Serpent and MCNP. Serpent and MCNP are in good agreement.

To investigate further, I then sample source particles only at the sphere centre (i.e. a point source). I tally the current passing through the sphere surface, and get the following results:

Point source:
OpenMC: 0.4613 (current in 1 MeV bin: 0.3344)
Serpent: 0.4648 (current in 1 MeV bin: 0.3345)
MCNP: 0.4653 (current in 1 MeV bin: 0.3348)

A much better agreement between all three codes is clearly seen. I attach a plot comparing the three spectra of the photons leaving the sphere (for both the point and the spherical source):

From the figure, the discrepancy in the 1-MeV peak (for the spherical source) is not really seen, but there are clearly some discrepancies between OpenMC and the other codes at lower energies (for both the point and the spherical source), perhaps due to differences in how the photon physics are handled? As can be seen, Serpent and MCNP are generally in good agreement.

What I cannot understand is why there is such a big difference in the tally results when using the spherical source distribution - a difference which is not observed when using the point source. This makes me think that there could be something going on with the sampling of the spherical source distribution in OpenMC.

Do you have any idea what could be the issue here? I’m using OpenMC version 0.13.0.

Thanks in advance!

An update: I attempted to define the spherical source as a custom source (explained here, but modified to a spherical source as explained here). Still using an Am241 sphere with a radius of 1 cm (which is also the source radius), I now get the following tallied current:

OpenMC: 0.6192793400000001
OpenMC counts in 1 MeV peak: 0.5011834100000001

which is obviously in better agreement with both Serpent and MCNP for the spherical source than when using “SphericalIndependent” to create the source. For now, I’ll stick with using the custom source method for my spherical source.

@mpres Sorry for the late follow-up on this. It turns out that to get a source uniformly sampled in a sphere, using a uniform distribution in theta is not the right thing to do and will result in points clustering toward the poles; it’s actually the cosine of theta that needs to be uniform. To be more precise, the distribution in theta itself would need to look like cos(theta) from theta=0 to π/2, but we don’t have a good way of representing that right now. Our development team right now is working on some changes to SphericalIndependent so that one can specify the distribution in cos(theta) instead of theta, which will make it easier to sample uniform points in a sphere. For now, the custom source route is the way to go so glad you were able to figure that out.

@mpres sorry for the late reply (I just did a pull request on the theta distribution.)

Another thing I think you should check is the distribution of the radius.

I think it should be:

r = openmc.stats.PowerLaw(0, sphere_radius, 3)

(3 instead of 1)

Otherwise you would have to many points in the center of the sphere.

That’s probably why the leaving current of OpenMC is smaller than the one of MCNP and Serpent.

Best

Christopher

Hi Christopher,

thanks for the reply. However, I think the sampling should be done in r^2 (as I have done it in the custom source), see for example here: Sampling in non-cartesian coordinate systems

Regards,
Markus

Hi Markus,

I am always getting very confused by the way the radius and theta needs to be chosen. I imagine the volume raises by the power of three so we also need to raise the number of sampled points by the power of three, to keep them uniform.

When reading your source (which is very good btw) it states that:

“Therefore, when using non-cartesian coordinates for Monte Carlo sampling, ascertain the Jacobian for the corresponding type of integral, change the variables of integration so that the corresponding integrand would be a constant, and sample those. So for example, in 3 dimensions, the sampling variables would be r3, cos(theta), and phi.”

So I still think r to the power of three should be distributed uniform, but please correct me if I’m wrong.

Best
Christopher

I, too, am always getting confused by this. The way I’ve come to understand it (which matches what @mpres suggested) is that you want P(r < R) for a fixed value of R to be proportional to r³. Since P(r < R) is essentially the CDF, that means the PDF (which is what we’re specifying in the SphericalIndependent class) needs to be proportional to r². Hence:

r_dist = openmc.stats.PowerLaw(0, sphere_radius, 2)

If you look at the implementation of the PowerLaw class in C++, you can see that for n=2, we take a uniform variable and then raise it to the power of 1/(n+1) = 1/3.

Yes, I am convinced. Thanks a lot :pray:

Sorry if I created confusion.