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)