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

Seem to be able to read the new hdf5 files and run at least an example problem (!!!!)

Browse source
This commit is contained in:
Adam Nelson 2016-08-21 13:43:07 -04:00
parent 63b3c6dbd0
commit 8fbd41eaf4
8 changed files with 61 additions and 78 deletions
Download patch file Download diff file Expand all files Collapse all files

4
examples/xml/basic/materials.xml
View file

@ -6,13 +6,15 @@
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U235" ao="1.0" />
<temperature>294K</temperature>
</material>
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H1" ao="2.0" />
<nuclide name="O16" ao="1.0" />
<sab name="c_H_in_H2O" xs="71t" />
<sab name="c_H_in_H2O"/>
<temperature>294K</temperature>
</material>
</materials>

9
openmc/data/neutron.py
View file

@ -399,7 +399,9 @@ class IncidentNeutron(EqualityMixin):
g.attrs['A'] = self.mass_number
g.attrs['metastable'] = self.metastable
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['kTs'] = self.kTs
ktg = g.create_group('kTs')
for i, temperature in enumerate(self.temperatures):
ktg.create_dataset(temperature, data=self.kTs[i])
# Write energy grid
eg = g.create_group('energy')
@ -457,7 +459,10 @@ class IncidentNeutron(EqualityMixin):
mass_number = group.attrs['A']
metastable = group.attrs['metastable']
atomic_weight_ratio = group.attrs['atomic_weight_ratio']
kTs = group.attrs['kTs'].tolist()
kTg = group['kTs']
kTs = []
for temp in kTg:
kTs.append(temp.value)
temperatures = [str(int(round(kT_to_K(kT)))) + "K" for kT in kTs]
data = cls(name, atomic_number, mass_number, metastable,

9
openmc/data/thermal.py
View file

@ -225,9 +225,11 @@ class ThermalScattering(EqualityMixin):
# Write basic data
g = f.create_group(self.name)
g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
g.attrs['kTs'] = self.kTs
g.attrs['zaids'] = self.zaids
g.attrs['secondary_mode'] = np.string_(self.secondary_mode)
ktg = g.create_group('kTs')
for i, temperature in enumerate(self.temperatures):
ktg.create_dataset(temperature, data=self.kTs[i])
for T in self.temperatures:
Tg = g.create_group(T)
@ -436,7 +438,10 @@ class ThermalScattering(EqualityMixin):
name = group.name[1:]
atomic_weight_ratio = group.attrs['atomic_weight_ratio']
kTs = group.attrs['kTs'].tolist()
kTg = group['kTs']
kTs = []
for temp in kTg:
kTs.append(temp.value)
temperatures = [str(int(round(kT_to_K(kT)))) + "K" for kT in kTs]
table = cls(name, atomic_weight_ratio, kTs)

3
src/global.F90
View file

@ -77,9 +77,6 @@ module global
! Dictionaries to look up cross sections and listings
type(DictCharInt) :: nuclide_dict
! Default xs identifier (e.g. 70c or 300K)
character(5):: default_xs
! ============================================================================
! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES

75
src/input_xml.F90
View file

@ -2159,16 +2159,9 @@ contains
&exist!")
end if
! Initialize default cross section variable
default_xs = ""
! Parse materials.xml file
call open_xmldoc(doc, filename)
! Copy default cross section if present
if (check_for_node(doc, "default_xs")) &
call get_node_value(doc, "default_xs", default_xs)
! Get pointer to list of XML <material>
call get_node_list(doc, "material", node_mat_list)
@ -2204,6 +2197,11 @@ contains
call get_node_value(node_mat, "name", mat % name)
end if
! Copy material temperature
if (check_for_node(node_mat, "temperature")) then
call get_node_value(node_mat, "temperature", mat % temperature)
end if
! =======================================================================
! READ AND PARSE <density> TAG
@ -2301,22 +2299,9 @@ contains
// trim(to_str(mat % id)))
end if
! Check for cross section
if (.not. check_for_node(node_nuc, "xs")) then
if (default_xs == '') then
call fatal_error("No cross section specified for macroscopic data &
& in material " // trim(to_str(mat % id)))
else
name = to_lower(trim(default_xs))
end if
end if
! store full name
call get_node_value(node_nuc, "name", temp_str)
if (check_for_node(node_nuc, "xs")) &
call get_node_value(node_nuc, "xs", name)
name = trim(temp_str) // "." // trim(name)
name = to_lower(name)
! store nuclide name
call get_node_value(node_nuc, "name", name)
name = trim(name)
! save name and density to list
call names % push_back(name)
@ -2345,16 +2330,6 @@ contains
// trim(to_str(mat % id)))
end if
! Check for cross section
if (.not. check_for_node(node_nuc, "xs")) then
if (default_xs == '') then
call fatal_error("No cross section specified for nuclide in &
&material " // trim(to_str(mat % id)))
else
name = to_lower(trim(default_xs))
end if
end if
! Check enforced isotropic lab scattering
if (run_CE) then
if (check_for_node(node_nuc, "scattering")) then
@ -2372,11 +2347,9 @@ contains
end if
end if
! store full name
call get_node_value(node_nuc, "name", temp_str)
if (check_for_node(node_nuc, "xs")) &
call get_node_value(node_nuc, "xs", name)
name = trim(temp_str) // "." // trim(name)
! store nuclide name
call get_node_value(node_nuc, "name", name)
name = trim(name)
! save name and density to list
call names % push_back(name)
@ -2424,18 +2397,6 @@ contains
end if
call get_node_value(node_ele, "name", name)
! Check for cross section
if (check_for_node(node_ele, "xs")) then
call get_node_value(node_ele, "xs", temp_str)
else
if (default_xs == '') then
call fatal_error("No cross section specified for nuclide in &
&material " // trim(to_str(mat % id)))
else
temp_str = to_lower(trim(default_xs))
end if
end if
! Check if no atom/weight percents were specified or if both atom and
! weight percents were specified
if (.not. check_for_node(node_ele, "ao") .and. &
@ -2581,14 +2542,11 @@ contains
call get_list_item(node_sab_list, j, node_sab)
! Determine name of S(a,b) table
if (.not. check_for_node(node_sab, "name") .or. &
.not. check_for_node(node_sab, "xs")) then
call fatal_error("Need to specify <name> and <xs> for S(a,b) &
&table.")
if (.not. check_for_node(node_sab, "name")) then
call fatal_error("Need to specify <name> for S(a,b) table.")
end if
call get_node_value(node_sab, "name", name)
call get_node_value(node_sab, "xs", temp_str)
name = trim(name) // "." // trim(temp_str)
name = trim(name)
mat % sab_names(j) = name
! Check that this nuclide is listed in the cross_sections.xml file
@ -5856,7 +5814,8 @@ contains
! Read nuclide data from HDF5
file_id = file_open(libraries(i_library) % path, 'r')
group_id = open_group(file_id, name)
call nuclides(i_nuclide) % from_hdf5(group_id)
call nuclides(i_nuclide) % from_hdf5(group_id, &
materials(i) % temperature)
call close_group(group_id)
call file_close(file_id)
@ -5988,7 +5947,7 @@ contains
! Read nuclide data from HDF5
file_id = file_open(libraries(i_library) % path, 'r')
group_id = open_group(file_id, name)
call resonant_nuc % from_hdf5(group_id)
call resonant_nuc % from_hdf5(group_id, '0K')
call close_group(group_id)
call file_close(file_id)

5
src/material_header.F90
View file

@ -8,11 +8,14 @@ module material_header
type Material
integer :: id ! unique identifier
character(len=104) :: name = "" ! User-defined name
character(len=104) :: name = "" ! User-defined name
integer :: n_nuclides ! number of nuclides
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
character(6) :: temperature ! Temperature of the material
! as presented in the HDF5 library;
! e.g., "300K"
! Energy grid information
integer :: n_grid ! # of union material grid points

21
src/nuclide_header.F90
View file

@ -180,9 +180,10 @@ module nuclide_header
end subroutine nuclide_clear
subroutine nuclide_from_hdf5(this, group_id)
subroutine nuclide_from_hdf5(this, group_id, temperature)
class(Nuclide), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
integer(HID_T), intent(in) :: group_id
character(6), intent(in) :: temperature
integer :: i
integer :: Z
@ -192,12 +193,14 @@ module nuclide_header
integer :: n_links
integer :: hdf5_err
integer(HID_T) :: urr_group, nu_group
integer(HID_T) :: energy_dset
integer(HID_T) :: energy_group, energy_dset
integer(HID_T) :: kT_group, kT_dset
integer(HID_T) :: rxs_group
integer(HID_T) :: rx_group
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
integer(HID_T) :: temp_dset
integer(SIZE_T) :: name_len, name_file_len
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
@ -217,15 +220,21 @@ module nuclide_header
call read_attribute(this % metastable, group_id, 'metastable')
this % zaid = 1000*Z + A + 400*this % metastable
call read_attribute(this % awr, group_id, 'atomic_weight_ratio')
call read_attribute(this % kT, group_id, 'temperature')
kT_group = open_group(group_id, 'kTs')
kT_dset = open_dataset(kT_group, temperature)
call read_dataset(this % kT, kT_dset)
call close_dataset(kT_dset)
call close_group(kT_group)
! Read energy grid
energy_dset = open_dataset(group_id, 'energy')
energy_group = open_group(group_id, 'energy')
energy_dset = open_dataset(energy_group, temperature)
call get_shape(energy_dset, dims)
this % n_grid = int(dims(1), 4)
allocate(this % energy(this % n_grid))
call read_dataset(this % energy, energy_dset)
call close_dataset(energy_dset)
call close_group(energy_group)
! Get MT values based on group names
rxs_group = open_group(group_id, 'reactions')
@ -243,7 +252,7 @@ module nuclide_header
do i = 1, size(this % reactions)
rx_group = open_group(rxs_group, 'reaction_' // trim(&
zero_padded(MTs % data(i), 3)))
call this % reactions(i) % from_hdf5(rx_group)
call this % reactions(i) % from_hdf5(rx_group, temperature)
call close_group(rx_group)
end do
call close_group(rxs_group)

13
src/reaction_header.F90
View file

@ -29,9 +29,10 @@ module reaction_header
contains
subroutine reaction_from_hdf5(this, group_id)
subroutine reaction_from_hdf5(this, group_id, temperature)
class(Reaction), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
character(6), intent(in) :: temperature
integer :: i
integer :: cm
@ -41,7 +42,7 @@ contains
integer :: n_links
integer :: hdf5_err
integer(HID_T) :: pgroup
integer(HID_T) :: xs
integer(HID_T) :: xs, xs_group
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: j
@ -49,16 +50,18 @@ contains
call read_attribute(this % Q_value, group_id, 'Q_value')
call read_attribute(this % MT, group_id, 'mt')
call read_attribute(this % threshold, group_id, 'threshold_idx')
call read_attribute(cm, group_id, 'center_of_mass')
this % scatter_in_cm = (cm == 1)
! Read cross section
xs = open_dataset(group_id, 'xs')
! Read cross section and threshold_idx data
xs_group = open_group(group_id, temperature)
call read_attribute(this % threshold, xs_group, 'threshold_idx')
xs = open_dataset(xs_group, 'xs')
call get_shape(xs, dims)
allocate(this % sigma(dims(1)))
call read_dataset(this % sigma, xs)
call close_dataset(xs)
call close_group(xs_group)
! Determine number of products
call h5gget_info_f(group_id, storage_type, n_links, max_corder, hdf5_err)