crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #777 from paulromano/thermal-bugs

Browse source 
Thermal scattering fixes and improvements
This commit is contained in:
Will Boyd 2017-01-18 08:04:49 -05:00 • committed by GitHub
commit 9eff5b8a2e
14 changed files with 165 additions and 98 deletions
Download patch file Download diff file Expand all files Collapse all files

2
.travis.yml
View file

@ -43,7 +43,7 @@ install: true
before_script:
- if [[ ! -e $HOME/nndc_hdf5/cross_sections.xml ]]; then
wget https://anl.box.com/shared/static/1ow880up90hgwynvh91twwfteameuwre.xz -O - | tar -C $HOME -xvJ;
wget https://anl.box.com/shared/static/oy85y1i3gboho82ifdtjlu1q4ozhd7el.xz -O - | tar -C $HOME -xvJ;
fi
- export OPENMC_CROSS_SECTIONS=$HOME/nndc_hdf5/cross_sections.xml

2
LICENSE
View file

@ -1,4 +1,4 @@
Copyright (c) 2011-2016 Massachusetts Institute of Technology
Copyright (c) 2011-2017 Massachusetts Institute of Technology
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

2
docs/source/conf.py
View file

@ -66,7 +66,7 @@ master_doc = 'index'
# General information about the project.
project = u'OpenMC'
copyright = u'2011-2016, Massachusetts Institute of Technology'
copyright = u'2011-2017, Massachusetts Institute of Technology'
# The version info for the project you're documenting, acts as replacement for
# |version| and |release|, also used in various other places throughout the

3
docs/source/devguide/styleguide.rst
View file

@ -153,6 +153,9 @@ Avoid extraneous whitespace in the following situations:
Yes: if (variable == 2) then
No: if ( variable==2 ) then
The structure component designator ``%`` should be surrounded by one space on
each side.
Do not leave trailing whitespace at the end of a line.
------

2
docs/source/license.rst
View file

@ -4,7 +4,7 @@
License Agreement
=================
Copyright © 2011-2016 Massachusetts Institute of Technology
Copyright © 2011-2017 Massachusetts Institute of Technology
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

2
man/man1/openmc.1
View file

@ -56,7 +56,7 @@ Indicates the default path to a directory containing windowed multipole data if
the user has not specified the <multipole_library> tag in
.I settings.xml\fP.
.SH LICENSE
Copyright \(co 2011-2016 Massachusetts Institute of Technology.
Copyright \(co 2011-2017 Massachusetts Institute of Technology.
.PP
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in

152
openmc/data/thermal.py
View file

@ -1,6 +1,7 @@
from collections import Iterable
from difflib import get_close_matches
from numbers import Real
import re
from warnings import warn
import numpy as np
@ -9,7 +10,7 @@ import h5py
import openmc.checkvalue as cv
from openmc.mixin import EqualityMixin
from . import HDF5_VERSION, HDF5_VERSION_MAJOR
from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV
from .data import K_BOLTZMANN, ATOMIC_SYMBOL, EV_PER_MEV, NATURAL_ABUNDANCE
from .ace import Table, get_table
from .angle_energy import AngleEnergy
from .function import Tabulated1D
@ -17,51 +18,80 @@ from .correlated import CorrelatedAngleEnergy
from openmc.stats import Discrete, Tabular
_THERMAL_NAMES = {'al': 'c_Al27', 'al27': 'c_Al27',
'be': 'c_Be',
'beo': 'c_BeO',
'bebeo': 'c_Be_in_BeO', 'be-o': 'c_Be_in_BeO', 'be/o': 'c_Be_in_BeO',
'benz': 'c_Benzine',
'cah': 'c_Ca_in_CaH2',
'dd2o': 'c_D_in_D2O', 'hwtr': 'c_D_in_D2O', 'hw': 'c_D_in_D2O',
'fe': 'c_Fe56', 'fe56': 'c_Fe56',
'graph': 'c_Graphite', 'grph': 'c_Graphite', 'gr': 'c_Graphite',
'hca': 'c_H_in_CaH2',
'hch2': 'c_H_in_CH2', 'poly': 'c_H_in_CH2', 'pol': 'c_H_in_CH2',
'hh2o': 'c_H_in_H2O', 'lwtr': 'c_H_in_H2O', 'lw': 'c_H_in_H2O',
'hzrh': 'c_H_in_ZrH', 'h-zr': 'c_H_in_ZrH', 'h/zr': 'c_H_in_ZrH', 'hzr': 'c_H_in_ZrH',
'lch4': 'c_liquid_CH4', 'lmeth': 'c_liquid_CH4',
'mg': 'c_Mg24',
'obeo': 'c_O_in_BeO', 'o-be': 'c_O_in_BeO', 'o/be': 'c_O_in_BeO',
'orthod': 'c_ortho_D', 'dortho': 'c_ortho_D',
'orthoh': 'c_ortho_H', 'hortho': 'c_ortho_H',
'ouo2': 'c_O_in_UO2', 'o2-u': 'c_O_in_UO2', 'o2/u': 'c_O_in_UO2',
'sio2': 'c_SiO2',
'parad': 'c_para_D', 'dpara': 'c_para_D',
'parah': 'c_para_H', 'hpara': 'c_para_H',
'sch4': 'c_solid_CH4', 'smeth': 'c_solid_CH4',
'uuo2': 'c_U_in_UO2', 'u-o2': 'c_U_in_UO2', 'u/o2': 'c_U_in_UO2',
'zrzrh': 'c_Zr_in_ZrH', 'zr-h': 'c_Zr_in_ZrH', 'zr/h': 'c_Zr_in_ZrH'}
_THERMAL_NAMES = {
'c_Al27': ('al', 'al27'),
'c_Be': ('be', 'be-metal'),
'c_BeO': ('beo',),
'c_Be_in_BeO': ('bebeo', 'be-o', 'be/o'),
'c_C6H6': ('benz', 'c6h6'),
'c_C_in_SiC': ('csic',),
'c_Ca_in_CaH2': ('cah',),
'c_D_in_D2O': ('dd2o', 'hwtr', 'hw'),
'c_Fe56': ('fe', 'fe56'),
'c_Graphite': ('graph', 'grph', 'gr'),
'c_H_in_CaH2': ('hcah2',),
'c_H_in_CH2': ('hch2', 'poly', 'pol'),
'c_H_in_CH4_liquid': ('lch4', 'lmeth'),
'c_H_in_CH4_solid': ('sch4', 'smeth'),
'c_H_in_H2O': ('hh2o', 'lwtr', 'lw'),
'c_H_in_H2O_solid': ('hice',),
'c_H_in_C5O2H8': ('lucite', 'c5o2h8'),
'c_H_in_YH2': ('hyh2',),
'c_H_in_ZrH': ('hzrh', 'h-zr', 'h/zr', 'hzr'),
'c_Mg24': ('mg', 'mg24'),
'c_O_in_BeO': ('obeo', 'o-be', 'o/be'),
'c_O_in_D2O': ('od2o',),
'c_O_in_H2O_ice': ('oice',),
'c_O_in_UO2': ('ouo2', 'o2-u', 'o2/u'),
'c_ortho_D': ('orthod', 'dortho'),
'c_ortho_H': ('orthoh', 'hortho'),
'c_Si_in_SiC': ('sisic',),
'c_SiO2_alpha': ('sio2', 'sio2a'),
'c_SiO2_beta': ('sio2b'),
'c_para_D': ('parad', 'dpara'),
'c_para_H': ('parah', 'hpara'),
'c_U_in_UO2': ('uuo2', 'u-o2', 'u/o2'),
'c_Y_in_YH2': ('yyh2',),
'c_Zr_in_ZrH': ('zrzrh', 'zr-h', 'zr/h')
}
def get_thermal_name(name):
"""Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O'"""
"""Get proper S(a,b) table name, e.g. 'HH2O' -> 'c_H_in_H2O'
if name in _THERMAL_NAMES.values():
Parameters
----------
name : str
Name of an ACE thermal scattering table
Returns
-------
str
GND-format thermal scattering name
"""
if name in _THERMAL_NAMES:
return name
elif name.lower() in _THERMAL_NAMES:
return _THERMAL_NAMES[name.lower()]
else:
# Make an educated guess?? This actually works well for
# JEFF-3.2 which stupidly uses names like lw00.32t,
# lw01.32t, etc. for different temperatures
matches = get_close_matches(
name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5)
if len(matches) > 0:
return _THERMAL_NAMES[matches[0]]
for proper_name, names in _THERMAL_NAMES.items():
if name.lower() in names:
return proper_name
else:
# OK, we give up. Just use the ACE name.
return 'c_' + name
# Make an educated guess?? This actually works well for
# JEFF-3.2 which stupidly uses names like lw00.32t,
# lw01.32t, etc. for different temperatures
for proper_name, names in _THERMAL_NAMES.items():
matches = get_close_matches(
name.lower(), names, cutoff=0.5)
if len(matches) > 0:
warn('Thermal scattering material "{}" is not recognized. '
'Assigning a name of {}.'.format(name, proper_name))
return proper_name
else:
# OK, we give up. Just use the ACE name.
warn('Thermal scattering material "{0}" is not recognized. '
'Assigning a name of c_{0}.'.format(name))
return 'c_' + name
class CoherentElastic(EqualityMixin):
@ -425,22 +455,7 @@ class ThermalScattering(EqualityMixin):
# Get new name that is GND-consistent
ace_name, xs = ace.name.split('.')
if name is None:
if ace_name.lower() in _THERMAL_NAMES:
name = _THERMAL_NAMES[ace_name.lower()]
else:
# Make an educated guess?? This actually works well for JEFF-3.2
# which stupidly uses names like lw00.32t, lw01.32t, etc. for
# different temperatures
matches = get_close_matches(
ace_name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5)
if len(matches) > 0:
name = _THERMAL_NAMES[matches[0]]
else:
# OK, we give up. Just use the ACE name.
name = 'c_' + ace.name
warn('Thermal scattering material "{}" is not recognized. '
'Assigning a name of {}.'.format(ace.name, name))
name = get_thermal_name(ace_name)
# Assign temperature to the running list
kTs = [ace.temperature*EV_PER_MEV]
@ -523,30 +538,47 @@ class ThermalScattering(EqualityMixin):
# Incoherent/coherent elastic scattering cross section
idx = ace.jxs[4]
n_mu = ace.nxs[6] + 1
if idx != 0:
n_energy = int(ace.xss[idx])
energy = ace.xss[idx + 1: idx + 1 + n_energy]*EV_PER_MEV
P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy]
if ace.nxs[5] == 4:
# Coherent elastic
table.elastic_xs[temperatures[0]] = CoherentElastic(
energy, P*EV_PER_MEV)
# Coherent elastic shouldn't have angular distributions listed
assert n_mu == 0
else:
# Incoherent elastic
table.elastic_xs[temperatures[0]] = Tabulated1D(energy, P)
# Angular distribution
n_mu = ace.nxs[6]
if n_mu != -1:
# Angular distribution
assert n_mu > 0
idx = ace.jxs[6]
table.elastic_mu_out[temperatures[0]] = \
ace.xss[idx:idx + n_energy * n_mu]
table.elastic_mu_out[temperatures[0]].shape = \
(n_energy, n_mu)
# Get relevant nuclides
# Get relevant nuclides -- NJOY only allows one to specify three
# nuclides that the S(a,b) table applies to. Thus, for all elements
# other than H and Fe, we automatically add all the naturally-occurring
# isotopes.
for zaid, awr in ace.pairs:
if zaid > 0:
Z, A = divmod(zaid, 1000)
table.nuclides.append(ATOMIC_SYMBOL[Z] + str(A))
element = ATOMIC_SYMBOL[Z]
if element in ['H', 'Fe']:
table.nuclides.append(element + str(A))
else:
if element + '0' not in table.nuclides:
table.nuclides.append(element + '0')
for isotope in sorted(NATURAL_ABUNDANCE):
if re.match(r'{}\d+'.format(element), isotope):
if isotope not in table.nuclides:
table.nuclides.append(isotope)
return table

20
scripts/openmc-get-nndc-data
View file

@ -114,15 +114,19 @@ for filename in glob.glob('nndc/293.6K/ENDF-B-VII.1-neutron-293.6K/*'):
shutil.move(filename, 'nndc/293.6K/')
# ==============================================================================
# EDIT GRAPHITE ZAID (6012 to 6000)
# FIX ZAID ASSIGNMENTS FOR VARIOUS S(A,B) TABLES
print('Changing graphite ZAID from 6012 to 6000')
graphite = os.path.join('nndc', 'tsl', 'graphite.acer')
with open(graphite) as fh:
text = fh.read()
text = text.replace('6012', '6000', 1)
with open(graphite, 'w') as fh:
fh.write(text)
def fix_zaid(table, old, new):
filename = os.path.join('nndc', 'tsl', table)
with open(filename, 'r') as fh:
text = fh.read()
text = text.replace(old, new, 1)
with open(filename, 'w') as fh:
fh.write(text)
print('Fixing ZAIDs for S(a,b) tables')
fix_zaid('bebeo.acer', '8016', ' 0')
fix_zaid('obeo.acer', '4009', ' 0')
# ==============================================================================
# PROMPT USER TO DELETE .TAR.GZ FILES

2
scripts/openmc-update-inputs
View file

@ -29,7 +29,7 @@ from shutil import move
import xml.etree.ElementTree as ET
import openmc.data
from openmc.data.thermal import _THERMAL_NAMES
description = "Update OpenMC's input XML files to the latest format."
epilog = """\

2
src/cross_section.F90
View file

@ -80,7 +80,7 @@ contains
j = j + 1
! Don't check for S(a,b) tables if there are no more left
if (j > mat % n_sab) check_sab = .false.
if (j > size(mat % i_sab_tables)) check_sab = .false.
end if
end if

30
src/input_xml.F90
View file

@ -1,5 +1,6 @@
module input_xml
use algorithm, only: find
use cmfd_input, only: configure_cmfd
use constants
use dict_header, only: DictIntInt, DictCharInt, ElemKeyValueCI
@ -2502,14 +2503,13 @@ contains
! Set number of S(a,b) tables
mat % n_sab = n_sab
! Allocate names and indices for nuclides and tables
! Allocate names and indices for nuclides and tables -- for now we
! allocate these as the number of S(a,b) tables listed. Since a single
! table might apply to multiple nuclides, they are resized later if a
! table is indeed applied to multiple nuclides.
allocate(mat % sab_names(n_sab))
allocate(mat % i_sab_nuclides(n_sab))
allocate(mat % i_sab_tables(n_sab))
! Initialize i_sab_nuclides
mat % i_sab_nuclides = NONE
do j = 1, n_sab
! Get pointer to S(a,b) table
call get_list_item(node_sab_list, j, node_sab)
@ -4990,6 +4990,8 @@ contains
integer :: m ! position for sorting
integer :: temp_nuclide ! temporary value for sorting
integer :: temp_table ! temporary value for sorting
type(VectorInt) :: i_sab_tables
type(VectorInt) :: i_sab_nuclides
do i = 1, size(materials)
! Skip materials with no S(a,b) tables
@ -5004,20 +5006,32 @@ contains
associate (sab => sab_tables(mat % i_sab_tables(k)))
FIND_NUCLIDE: do j = 1, size(mat % nuclide)
if (any(sab % nuclides == nuclides(mat % nuclide(j)) % name)) then
mat % i_sab_nuclides(k) = j
exit FIND_NUCLIDE
call i_sab_tables % push_back(k)
call i_sab_nuclides % push_back(j)
end if
end do FIND_NUCLIDE
end associate
! Check to make sure S(a,b) table matched a nuclide
if (mat % i_sab_nuclides(k) == NONE) then
if (find(i_sab_tables, k) == -1) then
call fatal_error("S(a,b) table " // trim(mat % &
sab_names(k)) // " did not match any nuclide on material " &
// trim(to_str(mat % id)))
end if
end do ASSIGN_SAB
! Update i_sab_tables and i_sab_nuclides
deallocate(mat % i_sab_tables)
m = i_sab_tables % size()
allocate(mat % i_sab_tables(m))
allocate(mat % i_sab_nuclides(m))
mat % i_sab_tables(:) = i_sab_tables % data(1:m)
mat % i_sab_nuclides(:) = i_sab_nuclides % data(1:m)
! Clear entries in vectors for next material
call i_sab_tables % clear()
call i_sab_nuclides % clear()
! If there are multiple S(a,b) tables, we need to make sure that the
! entries in i_sab_nuclides are sorted or else they won't be applied
! correctly in the cross_section module. The algorithm here is a simple

2
src/output.F90
View file

@ -67,7 +67,7 @@ contains
write(UNIT=OUTPUT_UNIT, FMT=*) &
' | The OpenMC Monte Carlo Code'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' Copyright | 2011-2016 Massachusetts Institute of Technology'
' Copyright | 2011-2017 Massachusetts Institute of Technology'
write(UNIT=OUTPUT_UNIT, FMT=*) &
' License | http://openmc.readthedocs.io/en/latest/license.html'
write(UNIT=OUTPUT_UNIT, FMT='(11X,"Version | ",I1,".",I1,".",I1)') &

40
src/physics.F90
View file

@ -546,6 +546,7 @@ contains
real(8) :: c_j, c_j1 ! cumulative probability
real(8) :: frac ! interpolation factor on outgoing energy
real(8) :: r1 ! RNG for outgoing energy
real(8) :: mu_left, mu_right ! adjacent mu values
i_temp = micro_xs(i_nuclide) % index_temp_sab
@ -737,23 +738,36 @@ contains
E = E_1 + (E - E_i1_1) * (E_J - E_1) / (E_i1_J - E_i1_1)
end if
! Find angular distribution for closest outgoing energy bin
if (r1 - c_j < c_j1 - r1) then
j = j
else
j = j + 1
end if
! Sample outgoing cosine bin
k = 1 + int(prn() * sab % n_inelastic_mu)
! Will use mu from the randomly chosen incoming and closest outgoing
! energy bins
mu = sab % inelastic_data(l) % mu(k, j)
! Rather than use the sampled discrete mu directly, it is smeared over
! a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
! discrete mu value itself.
associate (mu_l => sab % inelastic_data(l) % mu)
f = (r1 - c_j)/(c_j1 - c_j)
! Determine (k-1)th mu value
if (k == 1) then
mu_left = -ONE
else
mu_left = mu_l(k-1, j) + f*(mu_l(k-1, j+1) - mu_l(k-1,j))
end if
! Determine kth mu value
mu = mu_l(k, j) + f*(mu_l(k, j+1) - mu_l(k, j))
! Determine (k+1)th mu value
if (k == sab % n_inelastic_mu) then
mu_right = ONE - mu
else
mu_right = mu_l(k+1, j) + f*(mu_l(k+1, j+1) - mu_l(k+1,j)) - mu
end if
end associate
! Smear angle
mu = mu + min(mu - mu_left, mu_right - mu)*(prn() - HALF)
else
call fatal_error("Invalid secondary energy mode on S(a,b) table " &
// trim(sab_tables(i_sab) % name))
end if ! (inelastic secondary energy treatment)
end if ! (elastic or inelastic)
end associate

2
tests/test_salphabeta/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
8.331430E-01 3.074913E-03
8.465449E-01 5.743852E-03