Write constant values as strings in HDF5 formats. This includes things like the

run mode, tally estimator, filter type, tally scores, etc. Also do not write
nuclides in summary file.
This commit is contained in:
Paul Romano 2015-09-20 08:39:33 +07:00
parent 0f1cbdfd14
commit 9a0acdebc3
10 changed files with 232 additions and 521 deletions

View file

@ -1,129 +0,0 @@
"""Dictionaries of integer-to-string mappings from openmc/src/constants.F90"""
RUN_TYPES = {1: 'fixed source',
2: 'k-eigenvalue',
3: 'plot',
4: 'particle restart'}
SURFACE_TYPES = {1: 'x-plane',
2: 'y-plane',
3: 'z-plane',
4: 'plane',
5: 'x-cylinder',
6: 'y-cylinder',
7: 'z-cylinder',
8: 'sphere',
9: 'x-cone',
10: 'y-cone',
11: 'z-cone'}
BC_TYPES = {0: 'transmission',
1: 'vacuum',
2: 'reflective',
3: 'periodic'}
FILL_TYPES = {1: 'normal',
2: 'fill',
3: 'lattice'}
LATTICE_TYPES = {1: 'rectangular',
2: 'hexagonal'}
ESTIMATOR_TYPES = {1: 'analog',
2: 'tracklength',
3: 'collision'}
FILTER_TYPES = {1: 'universe',
2: 'material',
3: 'cell',
4: 'cellborn',
5: 'surface',
6: 'mesh',
7: 'energy',
8: 'energyout',
9: 'distribcell'}
SCORE_TYPES = {-1: 'flux',
-2: 'total',
-3: 'scatter',
-4: 'nu-scatter',
-5: 'scatter-n',
-6: 'scatter-pn',
-7: 'nu-scatter-n',
-8: 'nu-scatter-pn',
-9: 'transport',
-10: 'n1n',
-11: 'absorption',
-12: 'fission',
-13: 'nu-fission',
-14: 'kappa-fission',
-15: 'current',
-16: 'flux-yn',
-17: 'total-yn',
-18: 'scatter-yn',
-19: 'nu-scatter-yn',
-20: 'events',
1: '(n,total)',
2: '(n,elastic)',
4: '(n,level)',
11: '(n,2nd)',
16: '(n,2n)',
17: '(n,3n)',
18: '(n,fission)',
19: '(n,f)',
20: '(n,nf)',
21: '(n,2nf)',
22: '(n,na)',
23: '(n,n3a)',
24: '(n,2na)',
25: '(n,3na)',
28: '(n,np)',
29: '(n,n2a)',
30: '(n,2n2a)',
32: '(n,nd)',
33: '(n,nt)',
34: '(n,nHe-3)',
35: '(n,nd2a)',
36: '(n,nt2a)',
37: '(n,4n)',
38: '(n,3nf)',
41: '(n,2np)',
42: '(n,3np)',
44: '(n,n2p)',
45: '(n,npa)',
91: '(n,nc)',
101: '(n,disappear)',
102: '(n,gamma)',
103: '(n,p)',
104: '(n,d)',
105: '(n,t)',
106: '(n,3He)',
107: '(n,a)',
108: '(n,2a)',
109: '(n,3a)',
111: '(n,2p)',
112: '(n,pa)',
113: '(n,t2a)',
114: '(n,d2a)',
115: '(n,pd)',
116: '(n,pt)',
117: '(n,da)',
201: '(n,Xn)',
202: '(n,Xgamma)',
203: '(n,Xp)',
204: '(n,Xd)',
205: '(n,Xt)',
206: '(n,X3He)',
207: '(n,Xa)',
444: '(damage)',
649: '(n,pc)',
699: '(n,dc)',
749: '(n,tc)',
799: '(n,3Hec)',
849: '(n,tc)'}
SCORE_TYPES.update({MT: '(n,n' + str(MT-50) + ')' for MT in range(51,91)})
SCORE_TYPES.update({MT: '(n,p' + str(MT-600) + ')' for MT in range(600,649)})
SCORE_TYPES.update({MT: '(n,d' + str(MT-650) + ')' for MT in range(650,699)})
SCORE_TYPES.update({MT: '(n,t' + str(MT-700) + ')' for MT in range(700,749)})
SCORE_TYPES.update({MT: '(n,3He' + str(MT-750) + ')' for MT in range(750,649)})
SCORE_TYPES.update({MT: '(n,a' + str(MT-800) + ')' for MT in range(800,849)})

View file

@ -5,10 +5,12 @@ from numbers import Real, Integral
import numpy as np
from openmc import Mesh
from openmc.constants import *
from openmc.checkvalue import check_type, check_iterable_type, \
check_greater_than, _isinstance
_FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface',
'mesh', 'energy', 'energyout', 'distribcell']
class Filter(object):
"""A filter used to constrain a tally to a specific criterion, e.g. only tally
events when the particle is in a certain cell and energy range.
@ -109,7 +111,7 @@ class Filter(object):
def type(self, type):
if type is None:
self._type = type
elif type not in FILTER_TYPES.values():
elif type not in _FILTER_TYPES:
msg = 'Unable to set Filter type to "{0}" since it is not one ' \
'of the supported types'.format(type)
raise ValueError(msg)

View file

@ -1,7 +1,5 @@
import struct
from openmc.constants import RUN_TYPES
class Particle(object):
"""Information used to restart a specific particle that caused a simulation to
@ -74,7 +72,7 @@ class Particle(object):
@property
def run_mode(self):
return RUN_TYPES[self._f['run_mode'].value]
return self._f['run_mode'].value.decode()
@property
def uvw(self):

View file

@ -4,7 +4,6 @@ import sys
import numpy as np
import openmc
from openmc.constants import *
if sys.version > '3':
long = int
@ -317,7 +316,7 @@ class StatePoint(object):
@property
def run_mode(self):
return RUN_TYPES[self._f['run_mode'].value]
return self._f['run_mode'].value.decode()
@property
def seed(self):
@ -355,16 +354,14 @@ class StatePoint(object):
# Iterate over all Tallies
for tally_key in tally_keys:
# Read integer Tally estimator type code (analog, tracklength, or collision)
estimator_type = self._f['{0}{1}/estimator'.format(base, tally_key)].value
# Read the Tally size specifications
n_realizations = self._f['{0}{1}/n_realizations'.format(base, tally_key)].value
# Create Tally object and assign basic properties
tally = openmc.Tally(tally_key)
tally._statepoint = self
tally.estimator = ESTIMATOR_TYPES[estimator_type]
tally.estimator = self._f['{0}{1}/estimator'.format(
base, tally_key)].value.decode()
tally.num_realizations = n_realizations
# Read the number of Filters
@ -375,8 +372,8 @@ class StatePoint(object):
# Initialize all Filters
for j in range(1, n_filters+1):
# Read the integer Filter type code
filter_type = self._f['{0}{1}/type'.format(subbase, j)].value
# Read the Filter type
filter_type = self._f['{0}{1}/type'.format(subbase, j)].value.decode()
# Read the Filter offset
offset = self._f['{0}{1}/offset'.format(subbase, j)].value
@ -389,21 +386,21 @@ class StatePoint(object):
raise ValueError(msg)
# Read the bin values
if FILTER_TYPES[filter_type] in ['energy', 'energyout']:
if filter_type in ['energy', 'energyout']:
bins = self._f['{0}{1}/bins'.format(subbase, j)].value
elif FILTER_TYPES[filter_type] in ['mesh', 'distribcell']:
elif filter_type in ['mesh', 'distribcell']:
bins = self._f['{0}{1}/bins'.format(subbase, j)].value
else:
bins = self._f['{0}{1}/bins'.format(subbase, j)].value
# Create Filter object
filter = openmc.Filter(FILTER_TYPES[filter_type], bins)
filter = openmc.Filter(filter_type, bins)
filter.offset = offset
filter.num_bins = n_bins
if FILTER_TYPES[filter_type] == 'mesh':
if filter_type == 'mesh':
mesh_ids = self._f['tallies/meshes/ids'].value
mesh_keys = self._f['tallies/meshes/keys'].value
@ -427,9 +424,8 @@ class StatePoint(object):
tally.num_score_bins = n_score_bins
score_bins = self._f['{0}{1}/score_bins'.format(
scores = self._f['{0}{1}/scores'.format(
base, tally_key)].value
scores = [SCORE_TYPES[score] for score in score_bins]
n_user_scores = self._f['{0}{1}/n_user_score_bins'
.format(base, tally_key)].value
@ -447,6 +443,7 @@ class StatePoint(object):
# Add the scores to the Tally
for j, score in enumerate(scores):
score = score.decode()
# If this is a scattering moment, insert the scattering order
if '-n' in score:
score = score.replace('-n', '-' + moments[j].decode())

View file

@ -46,7 +46,6 @@ class Summary(object):
def _read_geometry(self):
# Read in and initialize the Materials and Geometry
self._read_nuclides()
self._read_materials()
self._read_surfaces()
self._read_cells()
@ -54,35 +53,6 @@ class Summary(object):
self._read_lattices()
self._finalize_geometry()
def _read_nuclides(self):
self.n_nuclides = self._f['nuclides/n_nuclides']
# Initialize dictionary for each Nuclide
# Keys - Nuclide ZAIDs
# Values - Nuclide objects
self.nuclides = {}
for key in self._f['nuclides'].keys():
if key == 'n_nuclides':
continue
index = self._f['nuclides'][key]['index'].value
alias = self._f['nuclides'][key]['alias'].value.decode()
zaid = self._f['nuclides'][key]['zaid'].value
# Read the Nuclide's name (e.g., 'H-1' or 'U-235')
name = alias.split('.')[0]
# Read the Nuclide's cross-section identifier (e.g., '70c')
xs = alias.split('.')[1]
# Initialize this Nuclide and add to global dictionary of Nuclides
if 'nat' in name:
self.nuclides[zaid] = openmc.Element(name=name, xs=xs)
else:
self.nuclides[zaid] = openmc.Nuclide(name=name, xs=xs)
self.nuclides[zaid].zaid = zaid
def _read_materials(self):
self.n_materials = self._f['n_materials'].value
@ -100,18 +70,14 @@ class Summary(object):
name = self._f['materials'][key]['name'].value.decode()
density = self._f['materials'][key]['atom_density'].value
nuc_densities = self._f['materials'][key]['nuclide_densities'][...]
nuclides = self._f['materials'][key]['nuclides'][...]
n_sab = self._f['materials'][key]['n_sab'].value
nuclides = self._f['materials'][key]['nuclides'].value
# Read the names of the S(a,b) tables for this Material
if n_sab > 0:
# Read the names of the S(a,b) tables for this Material and add them
if 'sab_names' in self._f['materials'][key]:
sab_tables = self._f['materials'][key]['sab_names'].value
sab_names = []
sab_xs = []
for sab_table in sab_tables:
name, xs = sab_table.decode().split('.')
sab_names.append(name)
sab_xs.append(xs)
material.add_s_alpha_beta(name, xs)
# Create the Material
material = openmc.Material(material_id=material_id, name=name)
@ -119,22 +85,21 @@ class Summary(object):
# Set the Material's density to g/cm3 - this is what is used in OpenMC
material.set_density(density=density, units='g/cm3')
# Add all Nuclides to the Material
for i, zaid in enumerate(nuclides):
nuclide = self.get_nuclide_by_zaid(zaid)
density = nuc_densities[i]
# Add all nuclides to the Material
for fullname, density in zip(nuclides, nuc_densities):
fullname = fullname.decode().strip()
name, xs = fullname.split('.')
if 'nat' in name:
nuclide = openmc.Element(name=name, xs=xs)
else:
nuclide = openmc.Nuclide(name=name, xs=xs)
if isinstance(nuclide, openmc.Nuclide):
material.add_nuclide(nuclide, percent=density, percent_type='ao')
elif isinstance(nuclide, openmc.Element):
material.add_element(nuclide, percent=density, percent_type='ao')
# Add S(a,b) table(s?) to the Material
for i in range(n_sab):
name = sab_names[i]
xs = sab_xs[i]
material.add_s_alpha_beta(name, xs)
# Add the Material to the global dictionary of all Materials
self.materials[index] = material
@ -540,9 +505,9 @@ class Summary(object):
tally = openmc.Tally(tally_id, tally_name)
# Read score metadata
score_bins = self._f['{0}/score_bins'.format(subbase)][...]
for score_bin in score_bins:
tally.add_score(openmc.SCORE_TYPES[score_bin])
scores = self._f['{0}/scores'.format(subbase)].value
for score in scores:
tally.add_score(score.decode())
num_score_bins = self._f['{0}/n_score_bins'.format(subbase)][...]
tally.num_score_bins = num_score_bins
@ -554,8 +519,7 @@ class Summary(object):
subsubbase = '{0}/filter {1}'.format(subbase, j)
# Read filter type (e.g., "cell", "energy", etc.)
filter_type_code = self._f['{0}/type'.format(subsubbase)].value
filter_type = openmc.FILTER_TYPES[filter_type_code]
filter_type = self._f['{0}/type'.format(subsubbase)].value.decode()
# Read the filter bins
num_bins = self._f['{0}/n_bins'.format(subsubbase)].value
@ -587,29 +551,6 @@ class Summary(object):
if self.opencg_geometry is None:
self.opencg_geometry = get_opencg_geometry(self.openmc_geometry)
def get_nuclide_by_zaid(self, zaid):
"""Return a Nuclide object given the 'zaid' identifier for the nuclide.
Parameters
----------
zaid : int
1000*Z + A, where Z is the atomic number of the nuclide and A is the
mass number. For example, the zaid for U-235 is 92235.
Returns
-------
nuclide : openmc.nuclide.Nuclide or None
Nuclide matching the specified zaid, or None if no matching object
is found.
"""
for index, nuclide in self.nuclides.items():
if nuclide._zaid == zaid:
return nuclide
return None
def get_material_by_id(self, material_id):
"""Return a Material object given the material id

View file

@ -4,7 +4,6 @@ from xml.etree import ElementTree as ET
import sys
from openmc.checkvalue import check_type, check_value, check_greater_than
from openmc.constants import BC_TYPES
if sys.version_info[0] >= 3:
basestring = str
@ -12,6 +11,8 @@ if sys.version_info[0] >= 3:
# A static variable for auto-generated Surface IDs
AUTO_SURFACE_ID = 10000
_BC_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic']
def reset_auto_surface_id():
global AUTO_SURFACE_ID
@ -106,7 +107,7 @@ class Surface(object):
@boundary_type.setter
def boundary_type(self, boundary_type):
check_type('boundary type', boundary_type, basestring)
check_value('boundary type', boundary_type, BC_TYPES.values())
check_value('boundary type', boundary_type, _BC_TYPES)
self._boundary_type = boundary_type
def __repr__(self):

View file

@ -71,6 +71,7 @@ contains
integer :: int_scalar
integer(HID_T) :: file_id
character(MAX_WORD_LEN) :: mode
! Write meessage
call write_message("Loading particle restart file " &
@ -86,7 +87,13 @@ contains
call read_dataset(file_id, 'gen_per_batch', gen_per_batch)
call read_dataset(file_id, 'current_gen', current_gen)
call read_dataset(file_id, 'n_particles', n_particles)
call read_dataset(file_id, 'run_mode', previous_run_mode)
call read_dataset(file_id, 'run_mode', mode)
select case (mode)
case ('k-eigenvalue')
previous_run_mode = MODE_EIGENVALUE
case ('fixed source')
previous_run_mode = MODE_FIXEDSOURCE
end select
call read_dataset(file_id, 'id', p%id)
call read_dataset(file_id, 'weight', p%wgt)
call read_dataset(file_id, 'energy', p%E)

View file

@ -46,7 +46,16 @@ contains
call write_dataset(file_id, 'gen_per_batch', gen_per_batch)
call write_dataset(file_id, 'current_gen', current_gen)
call write_dataset(file_id, 'n_particles', n_particles)
call write_dataset(file_id, 'run_mode', run_mode)
select case(run_mode)
case (MODE_FIXEDSOURCE)
call write_dataset(file_id, 'run_mode', 'fixed source')
case (MODE_EIGENVALUE)
call write_dataset(file_id, 'run_mode', 'k-eigenvalue')
case (MODE_PLOTTING)
call write_dataset(file_id, 'run_mode', 'plot')
case (MODE_PARTICLE)
call write_dataset(file_id, 'run_mode', 'particle restart')
end select
call write_dataset(file_id, 'id', p%id)
call write_dataset(file_id, 'weight', src%wgt)
call write_dataset(file_id, 'energy', src%E)

View file

@ -13,6 +13,7 @@ module state_point
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning
use global
use hdf5_interface
@ -50,6 +51,7 @@ contains
integer(HID_T) :: tallies_group, tally_group
integer(HID_T) :: meshes_group, mesh_group
integer(HID_T) :: filter_group
character(20), allocatable :: scores(:)
character(8), allocatable :: moment_names(:) ! names of moments (e.g, P3)
character(MAX_FILE_LEN) :: filename
type(StructuredMesh), pointer :: meshp
@ -90,7 +92,16 @@ contains
call write_dataset(file_id, "seed", seed)
! Write run information
call write_dataset(file_id, "run_mode", run_mode)
select case(run_mode)
case (MODE_FIXEDSOURCE)
call write_dataset(file_id, "run_mode", "fixed source")
case (MODE_EIGENVALUE)
call write_dataset(file_id, "run_mode", "k-eigenvalue")
case (MODE_PLOTTING)
call write_dataset(file_id, "run_mode", "plot")
case (MODE_PARTICLE)
call write_dataset(file_id, "run_mode", "particle restart")
end select
call write_dataset(file_id, "n_particles", n_particles)
call write_dataset(file_id, "n_batches", n_batches)
@ -214,7 +225,14 @@ contains
tally_group = create_group(tallies_group, "tally " // &
trim(to_str(tally%id)))
call write_dataset(tally_group, "estimator", tally%estimator)
select case(tally%estimator)
case (ESTIMATOR_ANALOG)
call write_dataset(tally_group, "estimator", "analog")
case (ESTIMATOR_TRACKLENGTH)
call write_dataset(tally_group, "estimator", "tracklength")
case (ESTIMATOR_COLLISION)
call write_dataset(tally_group, "estimator", "collision")
end select
call write_dataset(tally_group, "n_realizations", tally%n_realizations)
call write_dataset(tally_group, "n_filters", tally%n_filters)
@ -223,7 +241,28 @@ contains
filter_group = create_group(tally_group, "filter " // &
trim(to_str(j)))
call write_dataset(filter_group, "type", tally%filters(j)%type)
! Write name of type
select case (tally%filters(j)%type)
case(FILTER_UNIVERSE)
call write_dataset(filter_group, "type", "universe")
case(FILTER_MATERIAL)
call write_dataset(filter_group, "type", "material")
case(FILTER_CELL)
call write_dataset(filter_group, "type", "cell")
case(FILTER_CELLBORN)
call write_dataset(filter_group, "type", "cellborn")
case(FILTER_SURFACE)
call write_dataset(filter_group, "type", "surface")
case(FILTER_MESH)
call write_dataset(filter_group, "type", "mesh")
case(FILTER_ENERGYIN)
call write_dataset(filter_group, "type", "energy")
case(FILTER_ENERGYOUT)
call write_dataset(filter_group, "type", "energyout")
case(FILTER_DISTRIBCELL)
call write_dataset(filter_group, "type", "distribcell")
end select
call write_dataset(filter_group, "offset", tally%filters(j)%offset)
call write_dataset(filter_group, "n_bins", tally%filters(j)%n_bins)
if (tally%filters(j)%type == FILTER_ENERGYIN .or. &
@ -253,9 +292,59 @@ contains
deallocate(key_array)
call write_dataset(tally_group, "n_score_bins", tally%n_score_bins)
allocate(scores(size(tally%score_bins)))
do j = 1, size(tally%score_bins)
select case(tally%score_bins(j))
case (SCORE_FLUX)
scores(j) = "flux"
case (SCORE_TOTAL)
scores(j) = "total"
case (SCORE_SCATTER)
scores(j) = "scatter"
case (SCORE_NU_SCATTER)
scores(j) = "nu-scatter"
case (SCORE_SCATTER_N)
scores(j) = "scatter-n"
case (SCORE_SCATTER_PN)
scores(j) = "scatter-pn"
case (SCORE_NU_SCATTER_N)
scores(j) = "nu-scatter-n"
case (SCORE_NU_SCATTER_PN)
scores(j) = "nu-scatter-pn"
case (SCORE_TRANSPORT)
scores(j) = "transport"
case (SCORE_N_1N)
scores(j) = "n1n"
case (SCORE_ABSORPTION)
scores(j) = "absorption"
case (SCORE_FISSION)
scores(j) = "fission"
case (SCORE_NU_FISSION)
scores(j) = "nu-fission"
case (SCORE_KAPPA_FISSION)
scores(j) = "kappa-fission"
case (SCORE_CURRENT)
scores(j) = "current"
case (SCORE_FLUX_YN)
scores(j) = "flux-yn"
case (SCORE_TOTAL_YN)
scores(j) = "total-yn"
case (SCORE_SCATTER_YN)
scores(j) = "scatter-yn"
case (SCORE_NU_SCATTER_YN)
scores(j) = "nu-scatter-yn"
case (SCORE_EVENTS)
scores(j) = "events"
case default
scores(j) = reaction_name(tally%score_bins(j))
end select
end do
call write_dataset(tally_group, "scores", scores)
call write_dataset(tally_group, "score_bins", tally%score_bins)
call write_dataset(tally_group, "n_user_score_bins", tally%n_user_score_bins)
deallocate(scores)
! Write explicit moment order strings for each score bin
k = 1
allocate(moment_names(tally%n_score_bins))
@ -578,22 +667,15 @@ contains
subroutine load_state_point()
integer :: i, j
integer :: int_array(3)
integer :: curr_key
integer, allocatable :: id_array(:)
integer, allocatable :: key_array(:)
integer, allocatable :: temp_array(:)
integer :: i
integer :: int_array(3)
integer(HID_T) :: file_id
integer(HID_T) :: cmfd_group
integer(HID_T) :: tallies_group, tally_group
integer(HID_T) :: meshes_group, mesh_group
integer(HID_T) :: filter_group
real(8) :: real_array(3)
logical :: source_present
character(MAX_FILE_LEN) :: path_temp
character(19) :: current_time
type(StructuredMesh), pointer :: meshp
integer(HID_T) :: tallies_group
integer(HID_T) :: tally_group
real(8) :: real_array(3)
logical :: source_present
character(MAX_WORD_LEN) :: word
type(TallyObject), pointer :: tally
! Write message
@ -614,27 +696,17 @@ contains
&in OpenMC.")
end if
! Read OpenMC version
call read_dataset(file_id, "version_major", int_array(1))
call read_dataset(file_id, "version_minor", int_array(2))
call read_dataset(file_id, "version_release", int_array(3))
if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR &
.or. int_array(3) /= VERSION_RELEASE) then
if (master) call warning("State point file was created with a different &
&version of OpenMC.")
end if
! Read date and time
call read_dataset(file_id, "date_and_time", current_time)
! Read path to input
call read_dataset(file_id, "path", path_temp)
! Read and overwrite random number seed
call read_dataset(file_id, "seed", seed)
! Read and overwrite run information except number of batches
call read_dataset(file_id, "run_mode", run_mode)
call read_dataset(file_id, "run_mode", word)
select case(word)
case ('fixed source')
run_mode = MODE_FIXEDSOURCE
case ('k-eigenvalue')
run_mode = MODE_EIGENVALUE
end select
call read_dataset(file_id, "n_particles", n_particles)
call read_dataset(file_id, "n_batches", int_array(1))
@ -694,108 +766,6 @@ contains
end if
end if
! Read number of meshes
tallies_group = open_group(file_id, "tallies")
meshes_group = open_group(tallies_group, "meshes")
call read_dataset(meshes_group, "n_meshes", n_meshes)
if (n_meshes > 0) then
! Read list of mesh keys-> IDs
allocate(id_array(n_meshes))
allocate(key_array(n_meshes))
call read_dataset(meshes_group, "ids", id_array)
call read_dataset(meshes_group, "keys", key_array)
! Read and overwrite mesh information
MESH_LOOP: do i = 1, n_meshes
meshp => meshes(id_array(i))
curr_key = key_array(id_array(i))
mesh_group = open_group(meshes_group, "mesh " // &
trim(to_str(curr_key)))
call read_dataset(mesh_group, "id", meshp%id)
call read_dataset(mesh_group, "type", meshp%type)
call read_dataset(mesh_group, "n_dimension", meshp%n_dimension)
call read_dataset(mesh_group, "dimension", meshp%dimension)
call read_dataset(mesh_group, "lower_left", meshp%lower_left)
call read_dataset(mesh_group, "upper_right", meshp%upper_right)
call read_dataset(mesh_group, "width", meshp%width)
call close_group(mesh_group)
end do MESH_LOOP
deallocate(id_array)
deallocate(key_array)
end if
call close_group(meshes_group)
! Read and overwrite number of tallies
call read_dataset(tallies_group, "n_tallies", n_tallies)
! Read list of tally keys-> IDs
allocate(id_array(n_tallies))
allocate(key_array(n_tallies))
call read_dataset(tallies_group, "ids", id_array)
call read_dataset(tallies_group, "keys", key_array)
! Read in tally metadata
TALLY_METADATA: do i = 1, n_tallies
! Get pointer to tally
tally => tallies(i)
curr_key = key_array(id_array(i))
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(curr_key)))
call read_dataset(tally_group, "estimator", tally%estimator)
call read_dataset(tally_group, "n_realizations", tally%n_realizations)
call read_dataset(tally_group, "n_filters", tally%n_filters)
FILTER_LOOP: do j = 1, tally%n_filters
filter_group = open_group(tally_group, "filter " // trim(to_str(j)))
call read_dataset(filter_group, "type", tally%filters(j)%type)
call read_dataset(filter_group, "offset", tally%filters(j)%offset)
call read_dataset(filter_group, "n_bins", tally%filters(j)%n_bins)
if (tally%filters(j)%type == FILTER_ENERGYIN .or. &
tally%filters(j)%type == FILTER_ENERGYOUT) then
call read_dataset(filter_group, "bins", tally%filters(j)%real_bins)
else
call read_dataset(filter_group, "bins", tally%filters(j)%int_bins)
end if
call close_group(filter_group)
end do FILTER_LOOP
call read_dataset(tally_group, "n_nuclides", tally%n_nuclide_bins)
! Set up nuclide bin array and then read
allocate(temp_array(tally%n_nuclide_bins))
call read_dataset(tally_group, "nuclides", temp_array)
NUCLIDE_LOOP: do j = 1, tally%n_nuclide_bins
if (temp_array(j) > 0) then
tally%nuclide_bins(j) = temp_array(j)
else
tally%nuclide_bins(j) = temp_array(j)
end if
end do NUCLIDE_LOOP
deallocate(temp_array)
! Write number of score bins, score bins, user score bins
call read_dataset(tally_group, "n_score_bins", tally%n_score_bins)
call read_dataset(tally_group, "score_bins", tally%score_bins)
call read_dataset(tally_group, "n_user_score_bins", tally%n_user_score_bins)
call close_group(tally_group)
end do TALLY_METADATA
! Check to make sure source bank is present
if (path_source_point == path_state_point .and. .not. source_present) then
call fatal_error("Source bank must be contained in statepoint restart &
@ -808,13 +778,6 @@ contains
! Read number of realizations for global tallies
call read_dataset(file_id, "n_realizations", n_realizations, indep=.true.)
! Read number of global tallies
call read_dataset(file_id, "n_global_tallies", int_array(1), indep=.false.)
if (int_array(1) /= N_GLOBAL_TALLIES) then
call fatal_error("Number of global tallies does not match in state &
&point.")
end if
! Read global tally data
call read_dataset(file_id, "global_tallies", global_tallies)
@ -825,14 +788,12 @@ contains
! Read in sum and sum squared
if (int_array(1) == 1) then
TALLY_RESULTS: do i = 1, n_tallies
! Set pointer to tally
tally => tallies(i)
curr_key = key_array(id_array(i))
! Read sum and sum_sq for each bin
tally_group = open_group(tallies_group, "tally " // &
trim(to_str(curr_key)))
trim(to_str(tally%id)))
call read_dataset(tally_group, "results", tally%results)
call close_group(tally_group)
end do TALLY_RESULTS
@ -841,8 +802,6 @@ contains
call close_group(tallies_group)
end if
deallocate(id_array)
deallocate(key_array)
! Read source if in eigenvalue mode
if (run_mode == MODE_EIGENVALUE) then

View file

@ -62,7 +62,6 @@ contains
call write_geometry(file_id)
call write_materials(file_id)
call write_nuclides(file_id)
if (n_tallies > 0) then
call write_tallies(file_id)
end if
@ -238,16 +237,6 @@ contains
! Write coefficients for surface
call write_dataset(surface_group, "coefficients", s%coeffs)
! Write positive neighbors
if (allocated(s%neighbor_pos)) then
call write_dataset(surface_group, "neighbors_positive", s%neighbor_pos)
end if
! Write negative neighbors
if (allocated(s%neighbor_neg)) then
call write_dataset(surface_group, "neighbors_negative", s%neighbor_neg)
end if
! Write boundary condition
select case (s%bc)
case (BC_TRANSMIT)
@ -395,9 +384,10 @@ contains
subroutine write_materials(file_id)
integer(HID_T), intent(in) :: file_id
integer :: i
integer :: j
integer, allocatable :: zaids(:)
integer :: i
integer :: j
integer :: i_list
character(12), allocatable :: nucnames(:)
integer(HID_T) :: materials_group
integer(HID_T) :: material_group
type(Material), pointer :: m
@ -425,26 +415,22 @@ contains
"atom/b-cm")
! Copy ZAID for each nuclide to temporary array
allocate(zaids(m%n_nuclides))
allocate(nucnames(m%n_nuclides))
do j = 1, m%n_nuclides
zaids(j) = nuclides(m%nuclide(j))%zaid
i_list = nuclides(m%nuclide(j))%listing
nucnames(j) = xs_listings(i_list)%alias
end do
! Write temporary array to 'nuclides'
call write_dataset(material_group, "nuclides", zaids)
call write_dataset(material_group, "nuclides", nucnames)
! Deallocate temporary array
deallocate(zaids)
deallocate(nucnames)
! Write atom densities
call write_dataset(material_group, "nuclide_densities", m%atom_density)
! Write S(a,b) information if present
call write_dataset(material_group, "n_sab", m%n_sab)
if (m%n_sab > 0) then
call write_dataset(material_group, "i_sab_nuclides", m%i_sab_nuclides)
call write_dataset(material_group, "i_sab_tables", m%i_sab_tables)
call write_dataset(material_group, "sab_names", m%sab_names)
end if
@ -462,12 +448,13 @@ contains
subroutine write_tallies(file_id)
integer(HID_T), intent(in) :: file_id
integer :: i, j
integer :: i, j
integer, allocatable :: temp_array(:) ! nuclide bin array
integer(HID_T) :: tallies_group
integer(HID_T) :: mesh_group
integer(HID_T) :: tally_group
integer(HID_T) :: filter_group
character(20), allocatable :: scores(:)
type(StructuredMesh), pointer :: m
type(TallyObject), pointer :: t
@ -486,7 +473,6 @@ contains
! Write type and number of dimensions
call write_dataset(mesh_group, "type", m%type)
call write_dataset(mesh_group, "n_dimension", m%n_dimension)
! Write mesh information
call write_dataset(mesh_group, "dimension", m%dimension)
@ -521,9 +507,6 @@ contains
FILTER_LOOP: do j = 1, t%n_filters
filter_group = create_group(tally_group, "filter " // trim(to_str(j)))
! Write type of filter
call write_dataset(filter_group, "type", t%filters(j)%type)
! Write number of bins for this filter
call write_dataset(filter_group, "n_bins", t%filters(j)%n_bins)
@ -538,21 +521,23 @@ contains
! Write name of type
select case (t%filters(j)%type)
case(FILTER_UNIVERSE)
call write_dataset(filter_group, "type_name", "universe")
call write_dataset(filter_group, "type", "universe")
case(FILTER_MATERIAL)
call write_dataset(filter_group, "type_name", "material")
call write_dataset(filter_group, "type", "material")
case(FILTER_CELL)
call write_dataset(filter_group, "type_name", "cell")
call write_dataset(filter_group, "type", "cell")
case(FILTER_CELLBORN)
call write_dataset(filter_group, "type_name", "cellborn")
call write_dataset(filter_group, "type", "cellborn")
case(FILTER_SURFACE)
call write_dataset(filter_group, "type_name", "surface")
call write_dataset(filter_group, "type", "surface")
case(FILTER_MESH)
call write_dataset(filter_group, "type_name", "mesh")
call write_dataset(filter_group, "type", "mesh")
case(FILTER_ENERGYIN)
call write_dataset(filter_group, "type_name", "energy")
call write_dataset(filter_group, "type", "energy")
case(FILTER_ENERGYOUT)
call write_dataset(filter_group, "type_name", "energyout")
call write_dataset(filter_group, "type", "energyout")
case(FILTER_DISTRIBCELL)
call write_dataset(filter_group, "type", "distribcell")
end select
call close_group(filter_group)
@ -577,8 +562,58 @@ contains
! Write number of score bins
call write_dataset(tally_group, "n_score_bins", t%n_score_bins)
allocate(scores(size(t%score_bins)))
do j = 1, size(t%score_bins)
select case(t%score_bins(j))
case (SCORE_FLUX)
scores(j) = "flux"
case (SCORE_TOTAL)
scores(j) = "total"
case (SCORE_SCATTER)
scores(j) = "scatter"
case (SCORE_NU_SCATTER)
scores(j) = "nu-scatter"
case (SCORE_SCATTER_N)
scores(j) = "scatter-n"
case (SCORE_SCATTER_PN)
scores(j) = "scatter-pn"
case (SCORE_NU_SCATTER_N)
scores(j) = "nu-scatter-n"
case (SCORE_NU_SCATTER_PN)
scores(j) = "nu-scatter-pn"
case (SCORE_TRANSPORT)
scores(j) = "transport"
case (SCORE_N_1N)
scores(j) = "n1n"
case (SCORE_ABSORPTION)
scores(j) = "absorption"
case (SCORE_FISSION)
scores(j) = "fission"
case (SCORE_NU_FISSION)
scores(j) = "nu-fission"
case (SCORE_KAPPA_FISSION)
scores(j) = "kappa-fission"
case (SCORE_CURRENT)
scores(j) = "current"
case (SCORE_FLUX_YN)
scores(j) = "flux-yn"
case (SCORE_TOTAL_YN)
scores(j) = "total-yn"
case (SCORE_SCATTER_YN)
scores(j) = "scatter-yn"
case (SCORE_NU_SCATTER_YN)
scores(j) = "nu-scatter-yn"
case (SCORE_EVENTS)
scores(j) = "events"
case default
scores(j) = reaction_name(t%score_bins(j))
end select
end do
call write_dataset(tally_group, "scores", scores)
call write_dataset(tally_group, "score_bins", t%score_bins)
deallocate(scores)
call close_group(tally_group)
end do TALLY_METADATA
@ -586,115 +621,6 @@ contains
end subroutine write_tallies
!===============================================================================
! WRITE_NUCLIDES
!===============================================================================
subroutine write_nuclides(file_id)
integer(HID_T), intent(in) :: file_id
integer :: i, j
integer :: size_total
integer :: size_xs
integer :: size_angle
integer :: size_energy
integer(HID_T) :: nuclides_group, nuclide_group
integer(HID_T) :: reactions_group, rxn_group
type(Nuclide), pointer :: nuc
type(Reaction), pointer :: rxn
type(UrrData), pointer :: urr
nuclides_group = create_group(file_id, "nuclides")
! write number of nuclides
call write_dataset(nuclides_group, "n_nuclides", n_nuclides_total)
! Write information on each nuclide
NUCLIDE_LOOP: do i = 1, n_nuclides_total
nuc => nuclides(i)
nuclide_group = create_group(nuclides_group, nuc%name)
! Write internal OpenMC index for this nuclide
call write_dataset(nuclide_group, "index", i)
! Determine size of cross-sections
size_xs = (5 + nuc%n_reaction) * nuc%n_grid * 8
size_total = size_xs
! Write some basic attributes
call write_dataset(nuclide_group, "zaid", nuc%zaid)
call write_dataset(nuclide_group, "alias", xs_listings(nuc%listing)%alias)
call write_dataset(nuclide_group, "awr", nuc%awr)
call write_dataset(nuclide_group, "kT", nuc%kT)
call write_dataset(nuclide_group, "n_grid", nuc%n_grid)
call write_dataset(nuclide_group, "n_reactions", nuc%n_reaction)
call write_dataset(nuclide_group, "n_fission", nuc%n_fission)
call write_dataset(nuclide_group, "size_xs", size_xs)
! =======================================================================
! WRITE INFORMATION ON EACH REACTION
! Create overall group for reactions and close it
reactions_group = create_group(nuclide_group, "reactions")
RXN_LOOP: do j = 1, nuc%n_reaction
! Information on each reaction
rxn => nuc%reactions(j)
rxn_group = create_group(reactions_group, trim(reaction_name(rxn%MT)))
! Determine size of angle distribution
if (rxn%has_angle_dist) then
size_angle = rxn%adist%n_energy * 16 + size(rxn%adist%data) * 8
else
size_angle = 0
end if
! Determine size of energy distribution
if (rxn%has_energy_dist) then
size_energy = size(rxn%edist%data) * 8
else
size_energy = 0
end if
! Write information on reaction
call write_dataset(rxn_group, "Q_value", rxn%Q_value)
call write_dataset(rxn_group, "multiplicity", rxn%multiplicity)
call write_dataset(rxn_group, "threshold", rxn%threshold)
call write_dataset(rxn_group, "size_angle", size_angle)
call write_dataset(rxn_group, "size_energy", size_energy)
! Accumulate data size
size_total = size_total + size_angle + size_energy
call close_group(rxn_group)
end do RXN_LOOP
call close_group(reactions_group)
! =======================================================================
! WRITE INFORMATION ON URR PROBABILITY TABLES
if (nuc%urr_present) then
urr => nuc%urr_data
call write_dataset(nuclide_group, "urr_n_energy", urr%n_energy)
call write_dataset(nuclide_group, "urr_n_prob", urr%n_prob)
call write_dataset(nuclide_group, "urr_interp", urr%interp)
call write_dataset(nuclide_group, "urr_inelastic", urr%inelastic_flag)
call write_dataset(nuclide_group, "urr_absorption", urr%absorption_flag)
call write_dataset(nuclide_group, "urr_min_E", urr%energy(1))
call write_dataset(nuclide_group, "urr_max_E", urr%energy(urr%n_energy))
end if
! Write total memory used
call write_dataset(nuclide_group, "size_total", size_total)
call close_group(nuclide_group)
end do NUCLIDE_LOOP
call close_group(nuclides_group)
end subroutine write_nuclides
!===============================================================================
! WRITE_TIMING
!===============================================================================