Skip to content

Commit 291b45a

Browse files
nelsonagnelsonag
committed
Merge pull request #348 from paulromano/physics-and-nndc
Two physics improvements and NNDC data fixes
2 parents 2203084 + 8394e1e commit 291b45a

42 files changed

Lines changed: 6161 additions & 6103 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎data/cross_sections_nndc.xml‎

Lines changed: 20 additions & 20 deletions
Large diffs are not rendered by default.

‎data/get_nndc_data.py‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
import sys
88
import tarfile
99
import glob
10+
import hashlib
1011

1112
try:
1213
from urllib.request import urlopen
@@ -20,6 +21,8 @@
2021
baseUrl = 'http://www.nndc.bnl.gov/endf/b7.1/aceFiles/'
2122
files = ['ENDF-B-VII.1-neutron-293.6K.tar.gz',
2223
'ENDF-B-VII.1-tsl.tar.gz']
24+
checksums = ['9729a17eb62b75f285d8a7628ace1449',
25+
'e17d827c92940a30f22f096d910ea186']
2326
block_size = 16384
2427

2528
# ==============================================================================
@@ -65,6 +68,18 @@
6568
print('')
6669
filesComplete.append(f)
6770

71+
# ==============================================================================
72+
# VERIFY MD5 CHECKSUMS
73+
74+
print('Verifying MD5 checksums...')
75+
for f, checksum in zip(files, checksums):
76+
downloadsum = hashlib.md5(open(f, 'r').read()).hexdigest()
77+
if downloadsum != checksum:
78+
raise IOError("MD5 checksum for {} does not match. If this is your first "
79+
"time receiving this message, please re-run the script. "
80+
"Otherwise, please contact OpenMC developers by emailing "
81+
"openmc-users@googlegroups.com.")
82+
6883
# ==============================================================================
6984
# EXTRACT FILES FROM TGZ
7085

@@ -78,6 +93,10 @@
7893
print('Extracting {0}...'.format(f))
7994
tgz.extractall(path='nndc/' + suffix)
8095

96+
# Move ACE files down one level
97+
for filename in glob.glob('nndc/293.6K/ENDF-B-VII.1-neutron-293.6K/*'):
98+
shutil.move(filename, 'nndc/293.6K/')
99+
81100
#===============================================================================
82101
# EDIT GRAPHITE ZAID (6012 to 6000)
83102

‎src/CMakeLists.txt‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -347,7 +347,7 @@ foreach(test ${TESTS})
347347
elseif(${test} MATCHES "test_sourcepoint_restart")
348348
set(RESTART_FILE statepoint.07.h5 source.07.h5)
349349
elseif(${test} MATCHES "test_particle_restart_eigval")
350-
set(RESTART_FILE particle_12_842.h5)
350+
set(RESTART_FILE particle_12_616.h5)
351351
elseif(${test} MATCHES "test_particle_restart_fixed")
352352
set(RESTART_FILE particle_7_6144.h5)
353353
else(${test} MATCHES "test_statepoint_restart")
@@ -362,7 +362,7 @@ foreach(test ${TESTS})
362362
elseif(${test} MATCHES "test_sourcepoint_restart")
363363
set(RESTART_FILE statepoint.07.binary source.07.binary)
364364
elseif(${test} MATCHES "test_particle_restart_eigval")
365-
set(RESTART_FILE particle_12_842.binary)
365+
set(RESTART_FILE particle_12_616.binary)
366366
elseif(${test} MATCHES "test_particle_restart_fixed")
367367
set(RESTART_FILE particle_7_6144.binary)
368368
else(${test} MATCHES "test_statepoint_restart")

‎src/ace.F90‎

Lines changed: 34 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -704,7 +704,8 @@ subroutine read_reactions(nuc)
704704
integer :: LXS ! location of cross-section locators
705705
integer :: LOCA ! location of cross-section for given MT
706706
integer :: IE ! reaction's starting index on energy grid
707-
integer :: NE ! number of energies for reaction
707+
integer :: NE ! number of energies
708+
integer :: NR ! number of interpolation regions
708709
type(Reaction), pointer :: rxn => null()
709710
type(ListInt) :: MTs
710711

@@ -755,18 +756,40 @@ subroutine read_reactions(nuc)
755756
rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1)))
756757
rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0)
757758

758-
! If multiplicity is energy-dependent (absolute value > 100), set it based
759-
! on the MT value
759+
! Read energy-dependent multiplicities
760760
if (rxn % multiplicity > 100) then
761-
if (any(rxn%MT == [11, 16, 24, 30, 41])) then
762-
rxn % multiplicity = 2
763-
elseif (any(rxn%MT == [17, 25, 42])) then
764-
rxn % multiplicity = 3
765-
elseif (rxn%MT == 37) then
766-
rxn % multiplicity = 4
767-
else
768-
rxn % multiplicity = 1
761+
! Set flag and allocate space for Tab1 to store yield
762+
rxn % multiplicity_with_E = .true.
763+
allocate(rxn % multiplicity_E)
764+
765+
XSS_index = JXS(11) + rxn % multiplicity - 101
766+
NR = nint(XSS(XSS_index))
767+
rxn % multiplicity_E % n_regions = NR
768+
769+
! allocate space for ENDF interpolation parameters
770+
if (NR > 0) then
771+
allocate(rxn % multiplicity_E % nbt(NR))
772+
allocate(rxn % multiplicity_E % int(NR))
769773
end if
774+
775+
! read ENDF interpolation parameters
776+
XSS_index = XSS_index + 1
777+
if (NR > 0) then
778+
rxn % multiplicity_E % nbt = get_int(NR)
779+
rxn % multiplicity_E % int = get_int(NR)
780+
end if
781+
782+
! allocate space for yield data
783+
XSS_index = XSS_index + 2*NR
784+
NE = nint(XSS(XSS_index))
785+
rxn % multiplicity_E % n_pairs = NE
786+
allocate(rxn % multiplicity_E % x(NE))
787+
allocate(rxn % multiplicity_E % y(NE))
788+
789+
! read yield data
790+
XSS_index = XSS_index + 1
791+
rxn % multiplicity_E % x = get_real(NE)
792+
rxn % multiplicity_E % y = get_real(NE)
770793
end if
771794

772795
! read starting energy index

‎src/ace_header.F90‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -51,8 +51,10 @@ module ace_header
5151
integer :: MT ! ENDF MT value
5252
real(8) :: Q_value ! Reaction Q value
5353
integer :: multiplicity ! Number of secondary particles released
54+
type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield
5455
integer :: threshold ! Energy grid index of threshold
5556
logical :: scatter_in_cm ! scattering system in center-of-mass?
57+
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
5658
real(8), allocatable :: sigma(:) ! Cross section values
5759
logical :: has_angle_dist ! Angle distribution present?
5860
logical :: has_energy_dist ! Energy distribution present?

‎src/physics.F90‎

Lines changed: 26 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1291,6 +1291,7 @@ subroutine inelastic_scatter(nuc, rxn, E, uvw, mu, wgt)
12911291
real(8) :: E_in ! incoming energy
12921292
real(8) :: E_cm ! outgoing energy in center-of-mass
12931293
real(8) :: Q ! Q-value of reaction
1294+
real(8) :: yield ! neutron yield
12941295

12951296
! copy energy of neutron
12961297
E_in = E
@@ -1330,8 +1331,13 @@ subroutine inelastic_scatter(nuc, rxn, E, uvw, mu, wgt)
13301331
! change direction of particle
13311332
uvw = rotate_angle(uvw, mu)
13321333

1333-
! change weight of particle based on multiplicity
1334-
wgt = rxn % multiplicity * wgt
1334+
! change weight of particle based on yield
1335+
if (rxn % multiplicity_with_E) then
1336+
yield = interpolate_tab1(rxn % multiplicity_E, E_in)
1337+
else
1338+
yield = rxn % multiplicity
1339+
end if
1340+
wgt = yield * wgt
13351341

13361342
end subroutine inelastic_scatter
13371343

@@ -1572,6 +1578,7 @@ recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
15721578
real(8) :: E_max ! parameter for n-body dist
15731579
real(8) :: x, y, v ! intermediate variables for n-body dist
15741580
real(8) :: r1, r2, r3, r4, r5, r6
1581+
logical :: histogram_interp ! use histogram interpolation on incoming energy
15751582

15761583
! ==========================================================================
15771584
! SAMPLE ENERGY DISTRIBUTION IF THERE ARE MULTIPLE
@@ -1668,12 +1675,13 @@ recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
16681675
NR = int(edist % data(1))
16691676
NE = int(edist % data(2 + 2*NR))
16701677
if (NR == 1) then
1671-
if (master) call warning("Assuming linear-linear interpolation when &
1672-
&sampling continuous tabular distribution")
1678+
histogram_interp = (edist % data(3) == 1)
16731679
else if (NR > 1) then
16741680
! call write_particle_restart(p)
16751681
call fatal_error("Multiple interpolation regions not supported while &
16761682
&attempting to sample continuous tabular distribution.")
1683+
else
1684+
histogram_interp = .false.
16771685
end if
16781686

16791687
! find energy bin and calculate interpolation factor -- if the energy is
@@ -1693,11 +1701,15 @@ recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
16931701
end if
16941702

16951703
! Sample between the ith and (i+1)th bin
1696-
r2 = prn()
1697-
if (r > r2) then
1698-
l = i + 1
1699-
else
1704+
if (histogram_interp) then
17001705
l = i
1706+
else
1707+
r2 = prn()
1708+
if (r > r2) then
1709+
l = i + 1
1710+
else
1711+
l = i
1712+
end if
17011713
end if
17021714

17031715
! interpolation for energy E1 and EK
@@ -1776,10 +1788,12 @@ recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
17761788
end if
17771789

17781790
! Now interpolate between incident energy bins i and i + 1
1779-
if (l == i) then
1780-
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
1781-
else
1782-
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
1791+
if (.not. histogram_interp) then
1792+
if (l == i) then
1793+
E_out = E_1 + (E_out - E_i_1)*(E_K - E_1)/(E_i_K - E_i_1)
1794+
else
1795+
E_out = E_1 + (E_out - E_i1_1)*(E_K - E_1)/(E_i1_K - E_i1_1)
1796+
end if
17831797
end if
17841798

17851799
case (5)

‎tests/test_filter_cell/results_true.dat‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
k-combined:
2-
9.600453E-01 2.025740E-01
2+
1.093844E+00 1.626801E-02
33
tallies:
44
0.000000E+00
55
0.000000E+00
6-
1.161117E+01
7-
2.754585E+01
8-
2.434414E+00
9-
1.211549E+00
10-
3.447929E+01
11-
2.424749E+02
6+
1.517577E+01
7+
4.747271E+01
8+
3.151504E+00
9+
2.051857E+00
10+
4.536316E+01
11+
4.258781E+02

‎tests/test_filter_cellborn/results_true.dat‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,10 @@
11
k-combined:
2-
9.600453E-01 2.025740E-01
2+
1.093844E+00 1.626801E-02
33
tallies:
44
0.000000E+00
55
0.000000E+00
6-
5.474596E+01
7-
6.113816E+02
6+
7.449502E+01
7+
1.145793E+03
88
0.000000E+00
99
0.000000E+00
1010
0.000000E+00

‎tests/test_filter_energy/results_true.dat‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
k-combined:
2-
9.600453E-01 2.025740E-01
2+
1.093844E+00 1.626801E-02
33
tallies:
4-
2.878078E+01
5-
1.706092E+02
6-
4.244509E+01
7-
3.618897E+02
8-
5.296331E+01
9-
5.632082E+02
10-
9.415257E+00
11-
1.798132E+01
4+
2.687671E+01
5+
1.475192E+02
6+
4.148025E+01
7+
3.443331E+02
8+
5.223662E+01
9+
5.466362E+02
10+
1.050254E+01
11+
2.218108E+01

‎tests/test_filter_energyout/results_true.dat‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
k-combined:
2-
9.600453E-01 2.025740E-01
2+
1.093844E+00 1.626801E-02
33
tallies:
4-
2.904000E+01
5-
1.735606E+02
6-
4.173000E+01
7-
3.504547E+02
8-
5.212000E+01
9-
5.450908E+02
10-
6.720000E+00
11-
9.287600E+00
4+
2.740000E+01
5+
1.516786E+02
6+
4.163000E+01
7+
3.469463E+02
8+
5.041000E+01
9+
5.097599E+02
10+
6.990000E+00
11+
9.850700E+00

0 commit comments

Comments
 (0)