Continuous energy external source with constant probability for each energy bin

User Support
1

Hello, everyone. I am trying to simulate a fixed source problem with openmc. The external neutron source has a continuous energy which are separated into 47 groups. Each group has an identical neutron emission probability. According to the openmc manual, I chose openmc.stats.Tabular to set the energy distribution of the neutron source. Here is a part of codes for setting the external source:

# energy distribution
energy_bins =\
    [1.00000E-05, 1.00000E-01, 4.13990E-01, 8.76430E-01, 1.85540E+00, 5.04350E+00,
     1.06770E+01, 3.72660E+01, 1.01300E+02, 2.14450E+02, 4.54000E+02, 1.58460E+03,
     3.35460E+03, 7.10170E+03, 1.50340E+04, 2.18750E+04, 2.41760E+04, 2.60580E+04,
     3.18280E+04, 4.08680E+04, 6.73790E+04, 1.11090E+05, 1.83160E+05, 2.97210E+05,
     3.68830E+05, 4.97870E+05, 6.08100E+05, 7.42740E+05, 8.20850E+05, 1.00260E+06,
     1.35340E+06, 1.65300E+06, 1.92050E+06, 2.23130E+06, 2.34570E+06, 2.36530E+06,
     2.46600E+06, 2.72530E+06, 3.01190E+06, 3.67880E+06, 4.96590E+06, 6.06530E+06,
     7.40820E+06, 8.60710E+06, 1.00000E+07, 1.22140E+07, 1.41910E+07, 1.73320E+07]
# probabilities for each energy bin
energy_probs =\
    [1.44438E-11, 1.06857E-10, 2.52767E-10, 7.78361E-10, 4.01217E-09, 1.07794E-08,
     8.78665E-08, 3.61095E-07, 9.68270E-07, 2.98238E-06, 2.43405E-05, 5.97813E-05,
     1.83222E-04, 5.64378E-04, 6.27236E-04, 2.35381E-04, 2.00608E-04, 6.59333E-04,
     1.15149E-03, 4.06566E-03, 8.39881E-03, 1.71186E-02, 3.24986E-02, 2.27356E-02,
     4.34652E-02, 3.87843E-02, 4.80123E-02, 2.79514E-02, 6.39272E-02, 1.16353E-01,
     8.90702E-02, 7.06142E-02, 7.15504E-02, 2.36718E-02, 3.90518E-03, 1.93922E-02,
     4.53375E-02, 4.27965E-02, 7.51613E-02, 7.89060E-02, 2.92888E-02, 1.49788E-02,
     5.06871E-03, 2.17995E-03, 8.88027E-04, 1.39089E-04, 3.07600E-05]
   
# create a source for a fuel assembly
source = openmc.Source()
source.space = openmc.stats.Box((x_left, y_down, 0.0), (x_right, y_up, 190.0))
source.angle = openmc.stats.Isotropic()
source.energy = openmc.stats.Tabular(np.array(energy_bins), np.array(energy_probs), "histogram")                                      
source.strength = strength[count_assem]

However, I got an unreasonable result comparing to mcnp. I found that the external source sampled from openmc is different from that I input. Here are the figures of two neutron probability distributions:

image
image3000Ă—1800 215 KB

input_spectrum
input_spectrum3000Ă—1800 169 KB

I tried different interpolation schemes of openmc.stats.Tabular, but I cannot fix it. Any suggestions? Thank you very much.

2

Here are the codes to plot the energy probability density sampled from openmc:

# plot external source distribution
openmc.lib.init(output=False)
particles = openmc.lib.sample_external_source(n_samples=50000, prn_seed=1)
openmc.lib.finalize()

e_values = [particle.E for particle in particles]
e_array = np.array(e_values)

# count energy samples for each energy bin
counts, bin_edges = np.histogram(e_array, bins=energy_bins)

# plot probability density
plt.figure(figsize=(10, 6))
plt.bar(bin_edges[:-1], counts/50000, width=np.diff(bin_edges), edgecolor="black", align="edge")
plt.xscale("log")
plt.yscale("log")
plt.xlabel("Energy [eV]")
plt.ylabel("Probability")
plt.title("Energy spectrum from OpenMC samples")
plt.savefig("spectrum.jpg", dpi=300)
3

The thing to keep in mind is that the probabilities need to be given “per eV”, so if you haven’t already you should divide by the width of the energy bins. I suspect that is probably the issue here.

1 Like
4

Thanks for your reply. Dividing the probabbilities by energy bin width fixed this problem.