From 2622489b3bb5b10885fddad62ac12be1d8854c4f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 28 Jun 2017 12:50:01 -0500 Subject: [PATCH] Read photon cross sections at XML-read time, support particle type in source --- openmc/data/library.py | 17 +++++---- openmc/data/neutron.py | 1 + openmc/data/photon.py | 1 + openmc/data/thermal.py | 1 + openmc/source.py | 18 +++++++++- src/global.F90 | 7 ++++ src/input_xml.F90 | 76 +++++++++++++++++++++++++++++++++++++++- src/material_header.F90 | 1 + src/relaxng/settings.rnc | 1 + src/relaxng/settings.rng | 10 ++++++ src/source_header.F90 | 1 + 11 files changed, 125 insertions(+), 9 deletions(-) diff --git a/openmc/data/library.py b/openmc/data/library.py index c179f78f8c..2f33743a6e 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -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) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index 95f919369f..0894dac1e8 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -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 diff --git a/openmc/data/photon.py b/openmc/data/photon.py index e87a8792b4..d725c89751 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -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) diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 6e7b0fe02c..714c941f34 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -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 diff --git a/openmc/source.py b/openmc/source.py index 1aee435147..ca229638d3 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -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: diff --git a/src/global.F90 b/src/global.F90 index 2690ed1518..e42e04a8d7 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 597b4f4c0e..c0a91aecbd 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/material_header.F90 b/src/material_header.F90 index dc03734841..4132962407 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -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 diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 23a0cfd4b7..62caeed953 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -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 { diff --git a/src/relaxng/settings.rng b/src/relaxng/settings.rng index e75dc7c15e..c478bc5f24 100644 --- a/src/relaxng/settings.rng +++ b/src/relaxng/settings.rng @@ -230,6 +230,16 @@ + + + + + + + + + + diff --git a/src/source_header.F90 b/src/source_header.F90 index eefc007a47..763d8a7e95 100644 --- a/src/source_header.F90 +++ b/src/source_header.F90 @@ -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