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

Read photon cross sections at XML-read time, support particle type in source

Browse source
This commit is contained in:
Paul Romano 2017-06-28 12:50:01 -05:00
parent 645f076303
commit 2622489b3b
11 changed files with 125 additions and 9 deletions
Download patch file Download diff file Expand all files Collapse all files

17
openmc/data/library.py
View file

@ -48,14 +48,17 @@ class DataLibrary(EqualityMixin):
Path to the file to be registered.
"""
h5file = h5py.File(filename, 'r')
with h5py.File(filename, 'r') as h5file:
materials = []
filetype = 'neutron'
for name in h5file:
if name.startswith('c_'):
filetype = 'thermal'
materials.append(name)
materials = []
if 'filetype' in h5file.attrs:
filetype = h5file.attrs['filetype'].decode().lstrip('data_')
else:
filetype = 'neutron'
for name in h5file:
if name.startswith('c_'):
filetype = 'thermal'
materials.append(name)
library = {'path': filename, 'type': filetype, 'materials': materials}
self.libraries.append(library)

1
openmc/data/neutron.py
View file

@ -497,6 +497,7 @@ class IncidentNeutron(EqualityMixin):
# Open file and write version
f = h5py.File(path, mode, libver='latest')
f.attrs['filetype'] = np.string_('data_neutron')
f.attrs['version'] = np.array(HDF5_VERSION)
# Write basic data

1
openmc/data/photon.py
View file

@ -412,6 +412,7 @@ class IncidentPhoton(EqualityMixin):
"""
# Open file and write version
f = h5py.File(path, mode, libver='latest')
f.attrs['filetype'] = np.string_('data_photon')
if 'version' not in f.attrs:
f.attrs['version'] = np.array(HDF5_VERSION)

1
openmc/data/thermal.py
View file

@ -259,6 +259,7 @@ class ThermalScattering(EqualityMixin):
"""
# Open file and write version
f = h5py.File(path, mode, libver='latest')
f.attrs['filetype'] = np.string_('data_thermal')
f.attrs['version'] = np.array(HDF5_VERSION)
# Write basic data

18
openmc/source.py
View file

@ -24,6 +24,8 @@ class Source(object):
Source file from which sites should be sampled
strength : Real
Strength of the source
particle : {'neutron', 'photon'}
Source particle type
Attributes
----------
@ -37,10 +39,13 @@ class Source(object):
Source file from which sites should be sampled
strength : Real
Strength of the source
particle : {'neutron', 'photon'}
Source particle type
"""
def __init__(self, space=None, angle=None, energy=None, filename=None, strength=1.0):
def __init__(self, space=None, angle=None, energy=None, filename=None,
strength=1.0, particle='neutron'):
self._space = None
self._angle = None
self._energy = None
@ -55,6 +60,7 @@ class Source(object):
if filename is not None:
self.file = filename
self.strength = strength
self.particle = particle
@property
def file(self):
@ -76,6 +82,10 @@ class Source(object):
def strength(self):
return self._strength
@property
def particle(self):
return self._particle
@file.setter
def file(self, filename):
cv.check_type('source file', filename, string_types)
@ -102,6 +112,11 @@ class Source(object):
cv.check_greater_than('source strength', strength, 0.0, True)
self._strength = strength
@particle.setter
def particle(self, particle):
cv.check_value('source particle', particle, ['neutron', 'photon'])
self._particle = particle
def to_xml_element(self):
"""Return XML representation of the source
@ -113,6 +128,7 @@ class Source(object):
"""
element = ET.Element("source")
element.set("strength", str(self.strength))
element.set("particle", self.particle)
if self.file is not None:
element.set("file", self.file)
if self.space is not None:

7
src/global.F90
View file

@ -15,6 +15,7 @@ module global
use mesh_header, only: RegularMesh
use mgxs_header, only: Mgxs, MgxsContainer
use nuclide_header
use photon_header, only: PhotonInteraction
use plot_header, only: ObjectPlot
use sab_header, only: SAlphaBeta
use set_header, only: SetInt
@ -74,6 +75,7 @@ module global
! CROSS SECTION RELATED VARIABLES NEEDED REGARDLESS OF CE OR MG
integer :: n_nuclides_total ! Number of nuclide cross section tables
integer :: n_elements ! Number of photon cross section tables
! Cross section caches
type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide
@ -81,6 +83,7 @@ module global
! Dictionaries to look up cross sections and listings
type(DictCharInt) :: nuclide_dict
type(DictCharInt) :: element_dict
type(DictCharInt) :: library_dict
! Cross section libraries
@ -91,6 +94,7 @@ module global
! Cross section arrays
type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections
type(PhotonInteraction), allocatable :: elements(:) ! Photon cross sections
type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables
integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables
@ -115,6 +119,8 @@ module global
integer :: n_log_bins ! number of bins for logarithmic grid
real(8) :: log_spacing ! spacing on logarithmic grid
logical :: photon_transport = .false.
! ============================================================================
! MULTI-GROUP CROSS SECTION RELATED VARIABLES
@ -479,6 +485,7 @@ contains
end do
deallocate(nuclides)
end if
if (allocated(elements)) deallocate(elements)
if (allocated(res_scat_nuclides)) deallocate(res_scat_nuclides)

76
src/input_xml.F90
View file

@ -303,6 +303,7 @@ contains
if (n == 0) then
! Default source is isotropic point source at origin with Watt spectrum
allocate(external_source(1))
external_source % particle = NEUTRON
external_source % strength = ONE
allocate(SpatialPoint :: external_source(1) % space)
@ -335,6 +336,22 @@ contains
call get_node_value(node_source, "write_initial", write_initial_source)
end if
! Check for particle type
if (check_for_node(node_source, "particle")) then
call get_node_value(node_source, "particle", temp_str)
select case (to_lower(temp_str))
case ('neutron')
external_source(i) % particle = NEUTRON
case ('photon')
external_source(i) % particle = PHOTON
photon_transport = .true.
case default
call fatal_error('Unknown source particle type: ' // trim(temp_str))
end select
else
external_source(i) % particle = NEUTRON
end if
! Check for source strength
if (check_for_node(node_source, "strength")) then
call get_node_value(node_source, "strength", external_source(i)%strength)
@ -2208,9 +2225,11 @@ contains
integer :: n_sab ! number of sab tables for a material
integer :: i_library ! index in libraries array
integer :: index_nuclide ! index in nuclides
integer :: index_element ! index in elements
integer :: index_sab ! index in sab_tables
logical :: file_exists ! does materials.xml exist?
character(20) :: name ! name of nuclide, e.g. 92235.03c
character(20) :: name ! name of nuclide, e.g. U235
character(3) :: element ! name of element, e.g. Zr
character(MAX_WORD_LEN) :: units ! units on density
character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml
character(MAX_LINE_LEN) :: temp_str ! temporary string when reading
@ -2255,6 +2274,7 @@ contains
! Initialize count for number of nuclides/S(a,b) tables
index_nuclide = 0
index_element = 0
index_sab = 0
do i = 1, n_materials
@ -2481,6 +2501,7 @@ contains
mat % n_nuclides = n
allocate(mat % names(n))
allocate(mat % nuclide(n))
allocate(mat % element(n))
allocate(mat % atom_density(n))
allocate(mat % p0(n))
@ -2512,6 +2533,27 @@ contains
mat % nuclide(j) = nuclide_dict % get_key(to_lower(name))
end if
! If the corresponding element hasn't been encountered yet and photon
! transport will be used, we need to add its symbol to the element_dict
if (photon_transport) then
element = name(1:scan(name, '0123456789') - 1)
! Make sure photon cross section data is available
if (.not. library_dict % has_key(to_lower(element))) then
call fatal_error("Could not find element " // trim(element) &
// " in cross_sections data file!")
end if
if (.not. element_dict % has_key(element)) then
index_element = index_element + 1
mat % element(j) = index_element
call element_dict % add_key(element, index_element)
else
mat % element(j) = element_dict % get_key(element)
end if
end if
! Copy name and atom/weight percent
mat % names(j) = name
mat % atom_density(j) = densities % data(j)
@ -2610,6 +2652,7 @@ contains
! Set total number of nuclides and S(a,b) tables
n_nuclides_total = index_nuclide
n_elements = index_element
n_sab_tables = index_sab
! Close materials XML file
@ -4854,6 +4897,8 @@ contains
libraries(i) % type = LIBRARY_NEUTRON
case ('thermal')
libraries(i) % type = LIBRARY_THERMAL
case ('photon')
libraries(i) % type = LIBRARY_PHOTON
end select
else
call fatal_error("Missing library type")
@ -5230,14 +5275,18 @@ contains
integer :: i, j
integer :: i_library
integer :: i_nuclide
integer :: i_element
integer :: i_sab
integer(HID_T) :: file_id
integer(HID_T) :: group_id
logical :: mp_found ! if windowed multipole libraries were found
character(MAX_WORD_LEN) :: name
character(3) :: element
type(SetChar) :: already_read
type(SetChar) :: element_already_read
allocate(nuclides(n_nuclides_total))
allocate(elements(n_elements))
allocate(sab_tables(n_sab_tables))
! Read cross sections
@ -5279,6 +5328,31 @@ contains
! Add name and alias to dictionary
call already_read % add(name)
! Check if elemental data has been read, if needed
element = name(1:scan(name, '0123456789') - 1)
if (photon_transport) then
if (.not. element_already_read % contains(element)) then
! Read photon interaction data from HDF5 photon library
i_library = library_dict % get_key(to_lower(element))
i_element = element_dict % get_key(element)
call write_message('Reading ' // trim(element) // ' from ' // &
trim(libraries(i_library) % path), 6)
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
call check_data_version(file_id)
! Read element data from HDF5
group_id = open_group(file_id, element)
call elements(i_element) % from_hdf5(group_id)
call close_group(group_id)
call file_close(file_id)
! Add element to set
call element_already_read % add(element)
end if
end if
! Read multipole file into the appropriate entry on the nuclides array
if (temperature_multipole) call read_multipole_data(i_nuclide)
end if

1
src/material_header.F90
View file

@ -14,6 +14,7 @@ module material_header
character(len=104) :: name = "" ! User-defined name
integer :: n_nuclides ! number of nuclides
integer, allocatable :: nuclide(:) ! index in nuclides array
integer, allocatable :: element(:) ! index in elements array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
real(8) :: density_gpcc ! total density in g/cm^3

1
src/relaxng/settings.rnc
View file

@ -55,6 +55,7 @@ element settings {
element source {
grammar {
start =
(element particle { xsd:string } | attribute particle { xsd:string })? &
(element strength { xsd:double } | attribute strength { xsd:double })? &
(element file { xsd:string } | attribute file { xsd:string })? &
element space {

10
src/relaxng/settings.rng
View file

@ -230,6 +230,16 @@
<grammar>
<start>
<interleave>
<optional>
<choice>
<element name="particle">
<data type="string"/>
</element>
<attribute name="particle">
<data type="string"/>
</attribute>
</choice>
</optional>
<optional>
<choice>
<element name="strength">

1
src/source_header.F90
View file

@ -11,6 +11,7 @@ module source_header
!===============================================================================
type SourceDistribution
integer :: particle ! particle type
real(8) :: strength ! source strength
class(SpatialDistribution), allocatable :: space ! spatial distribution
class(UnitSphereDistribution), allocatable :: angle ! angle distribution