mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-29 06:35:48 -04:00
Fixed merge conflicts with develop branch
This commit is contained in:
Will Boyd
commit
e05ab41527
99 changed files with 5737 additions and 3263 deletions
|
|
@ -1,123 +0,0 @@
|
|||
"""Dictionaries of integer-to-string mappings from openmc/src/constants.F90"""
|
||||
|
||||
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'}
|
||||
|
||||
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)})
|
||||
|
|
@ -15,6 +15,9 @@ if sys.version_info[0] >= 3:
|
|||
basestring = str
|
||||
|
||||
|
||||
_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.
|
||||
|
|
@ -135,7 +138,9 @@ class Filter(object):
|
|||
|
||||
@type.setter
|
||||
def type(self, type):
|
||||
if type not in FILTER_TYPES.values():
|
||||
if type is None:
|
||||
self._type = type
|
||||
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)
|
||||
|
|
|
|||
|
|
@ -56,7 +56,7 @@ class Mesh(object):
|
|||
# Initialize Mesh class attributes
|
||||
self.id = mesh_id
|
||||
self.name = name
|
||||
self._type = 'rectangular'
|
||||
self._type = 'regular'
|
||||
self._dimension = None
|
||||
self._lower_left = None
|
||||
self._upper_right = None
|
||||
|
|
@ -158,7 +158,7 @@ class Mesh(object):
|
|||
cv.check_type('type for mesh ID="{0}"'.format(self._id),
|
||||
meshtype, basestring)
|
||||
cv.check_value('type for mesh ID="{0}"'.format(self._id),
|
||||
meshtype, ['rectangular', 'hexagonal'])
|
||||
meshtype, ['regular'])
|
||||
self._type = meshtype
|
||||
|
||||
@dimension.setter
|
||||
|
|
|
|||
|
|
@ -246,9 +246,6 @@ class MultiGroupXS(object):
|
|||
|
||||
cv.check_type('statepoint', statepoint, openmc.statepoint.StatePoint)
|
||||
|
||||
# Ensure that tally metadata has been loaded from the statepoint file
|
||||
statepoint.read_results()
|
||||
|
||||
# Create Tallies to search for in StatePoint
|
||||
self.create_tallies()
|
||||
|
||||
|
|
@ -1304,4 +1301,4 @@ class Chi(MultiGroupXS):
|
|||
nu_fission_out = self.tallies['nu-fission-out']
|
||||
self._xs_tally = nu_fission_out / nu_fission_in
|
||||
self._xs_tally._mean = np.nan_to_num(self.xs_tally.mean)
|
||||
self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev)
|
||||
self._xs_tally._std_dev = np.nan_to_num(self.xs_tally.std_dev)
|
||||
|
|
|
|||
|
|
@ -12,10 +12,6 @@ class Particle(object):
|
|||
|
||||
Attributes
|
||||
----------
|
||||
filetype : int
|
||||
Integer indicating the file type
|
||||
revision : int
|
||||
Revision of the particle restart format
|
||||
current_batch : int
|
||||
The batch containing the particle
|
||||
gen_per_batch : int
|
||||
|
|
@ -43,28 +39,52 @@ class Particle(object):
|
|||
import h5py
|
||||
self._f = h5py.File(filename, 'r')
|
||||
|
||||
# Read all metadata
|
||||
self._read_data()
|
||||
# Ensure filetype and revision are correct
|
||||
if 'filetype' not in self._f or self._f[
|
||||
'filetype'].value.decode() != 'particle restart':
|
||||
raise IOError('{} is not a particle restart file.'.format(filename))
|
||||
if self._f['revision'].value != 1:
|
||||
raise IOError('Particle restart file has a file revision of {} '
|
||||
'which is not consistent with the revision this '
|
||||
'version of OpenMC expects ({}).'.format(
|
||||
self._f['revision'].value, 1))
|
||||
|
||||
def _read_data(self):
|
||||
# Read filetype
|
||||
self.filetype = self._f['filetype'].value
|
||||
@property
|
||||
def current_batch(self):
|
||||
return self._f['current_batch'].value
|
||||
|
||||
# Read statepoint revision
|
||||
self.revision = self._f['revision'].value
|
||||
@property
|
||||
def current_gen(self):
|
||||
return self._f['current_gen'].value
|
||||
|
||||
# Read current batch
|
||||
self.current_batch = self._f['current_batch'].value
|
||||
@property
|
||||
def energy(self):
|
||||
return self._f['energy'].value
|
||||
|
||||
# Read run information
|
||||
self.gen_per_batch = self._f['gen_per_batch'].value
|
||||
self.current_gen = self._f['current_gen'].value
|
||||
self.n_particles = self._f['n_particles'].value
|
||||
self.run_mode = self._f['run_mode'].value
|
||||
@property
|
||||
def gen_per_batch(self):
|
||||
return self._f['gen_per_batch'].value
|
||||
|
||||
# Read particle properties
|
||||
self.id = self._f['id'].value
|
||||
self.weight = self._f['weight'].value
|
||||
self.energy = self._f['energy'].value
|
||||
self.xyz = self._f['xyz'].value
|
||||
self.uvw = self._f['uvw'].value
|
||||
@property
|
||||
def id(self):
|
||||
return self._f['id'].value
|
||||
|
||||
@property
|
||||
def n_particles(self):
|
||||
return self._f['n_particles'].value
|
||||
|
||||
@property
|
||||
def run_mode(self):
|
||||
return self._f['run_mode'].value.decode()
|
||||
|
||||
@property
|
||||
def uvw(self):
|
||||
return self._f['uvw'].value
|
||||
|
||||
@property
|
||||
def weight(self):
|
||||
return self._f['weight'].value
|
||||
|
||||
@property
|
||||
def xyz(self):
|
||||
return self._f['xyz'].value
|
||||
|
|
|
|||
|
|
@ -2,62 +2,13 @@ import copy
|
|||
import sys
|
||||
|
||||
import numpy as np
|
||||
import scipy.stats
|
||||
|
||||
import openmc
|
||||
from openmc.constants import *
|
||||
|
||||
if sys.version > '3':
|
||||
long = int
|
||||
|
||||
|
||||
class SourceSite(object):
|
||||
"""A single source site produced from fission.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
weight : float
|
||||
Weight of the particle arising from the site
|
||||
xyz : list of float
|
||||
Cartesian coordinates of the site
|
||||
uvw : list of float
|
||||
Directional cosines for particles emerging from the site
|
||||
E : float
|
||||
Energy of the emerging particle in MeV
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self):
|
||||
self._weight = None
|
||||
self._xyz = None
|
||||
self._uvw = None
|
||||
self._E = None
|
||||
|
||||
def __repr__(self):
|
||||
string = 'SourceSite\n'
|
||||
string += '{0: <16}{1}{2}\n'.format('\tweight', '=\t', self._weight)
|
||||
string += '{0: <16}{1}{2}\n'.format('\tE', '=\t', self._E)
|
||||
string += '{0: <16}{1}{2}\n'.format('\t(x,y,z)', '=\t', self._xyz)
|
||||
string += '{0: <16}{1}{2}\n'.format('\t(u,v,w)', '=\t', self._uvw)
|
||||
return string
|
||||
|
||||
@property
|
||||
def weight(self):
|
||||
return self._weight
|
||||
|
||||
@property
|
||||
def xyz(self):
|
||||
return self._xyz
|
||||
|
||||
@property
|
||||
def uvw(self):
|
||||
return self._uvw
|
||||
|
||||
@property
|
||||
def E(self):
|
||||
return self._E
|
||||
|
||||
|
||||
class StatePoint(object):
|
||||
"""State information on a simulation at a certain point in time (at the end of a
|
||||
given batch). Statepoints can be used to analyze tally results as well as
|
||||
|
|
@ -65,28 +16,74 @@ class StatePoint(object):
|
|||
|
||||
Attributes
|
||||
----------
|
||||
cmfd_on : bool
|
||||
Indicate whether CMFD is active
|
||||
cmfd_balance : ndarray
|
||||
Residual neutron balance for each batch
|
||||
cmfd_dominance
|
||||
Dominance ratio for each batch
|
||||
cmfd_entropy : ndarray
|
||||
Shannon entropy of CMFD fission source for each batch
|
||||
cmfd_indices : ndarray
|
||||
Number of CMFD mesh cells and energy groups. The first three indices
|
||||
correspond to the x-, y-, and z- spatial directions and the fourth index
|
||||
is the number of energy groups.
|
||||
cmfd_srccmp : ndarray
|
||||
Root-mean-square difference between OpenMC and CMFD fission source for
|
||||
each batch
|
||||
cmfd_src : ndarray
|
||||
CMFD fission source distribution over all mesh cells and energy groups.
|
||||
current_batch : int
|
||||
Number of batches simulated
|
||||
date_and_time : str
|
||||
Date and time when simulation began
|
||||
entropy : ndarray
|
||||
Shannon entropy of fission source at each batch
|
||||
gen_per_batch : int
|
||||
Number of fission generations per batch
|
||||
global_tallies : ndarray of compound datatype
|
||||
Global tallies for k-effective estimates and leakage. The compound
|
||||
datatype has fields 'name', 'sum', 'sum_sq', 'mean', and 'std_dev'.
|
||||
k_combined : list
|
||||
Combined estimator for k-effective and its uncertainty
|
||||
n_particles : int
|
||||
Number of particles per generation
|
||||
k_col_abs : float
|
||||
Cross-product of collision and absorption estimates of k-effective
|
||||
k_col_tra : float
|
||||
Cross-product of collision and tracklength estimates of k-effective
|
||||
k_abs_tra : float
|
||||
Cross-product of absorption and tracklength estimates of k-effective
|
||||
k_generation : ndarray
|
||||
Estimate of k-effective for each batch/generation
|
||||
meshes : dict
|
||||
Dictionary whose keys are mesh IDs and whose values are Mesh objects
|
||||
n_batches : int
|
||||
Number of batches
|
||||
current_batch :
|
||||
Number of batches simulated
|
||||
results : bool
|
||||
Indicate whether tally results have been read
|
||||
source : ndarray of SourceSite
|
||||
Array of source sites
|
||||
with_summary : bool
|
||||
Indicate whether statepoint data has been linked against a summary file
|
||||
n_inactive : int
|
||||
Number of inactive batches
|
||||
n_particles : int
|
||||
Number of particles per generation
|
||||
n_realizations : int
|
||||
Number of tally realizations
|
||||
path : str
|
||||
Working directory for simulation
|
||||
run_mode : str
|
||||
Simulation run mode, e.g. 'k-eigenvalue'
|
||||
seed : int
|
||||
Pseudorandom number generator seed
|
||||
source : ndarray of compound datatype
|
||||
Array of source sites. The compound datatype has fields 'wgt', 'xyz',
|
||||
'uvw', and 'E' corresponding to the weight, position, direction, and
|
||||
energy of the source site.
|
||||
source_present : bool
|
||||
Indicate whether source sites are present
|
||||
tallies : dict
|
||||
Dictionary whose keys are tally IDs and whose values are Tally objects
|
||||
tallies_present : bool
|
||||
Indicate whether user-defined tallies are present
|
||||
global_tallies : ndarray
|
||||
Global tallies and their uncertainties
|
||||
n_realizations : int
|
||||
Number of tally realizations
|
||||
version: tuple of int
|
||||
Version of OpenMC
|
||||
with_summary : bool
|
||||
Indicate whether statepoint data has been linked against a summary file
|
||||
|
||||
"""
|
||||
|
||||
|
|
@ -94,474 +91,375 @@ class StatePoint(object):
|
|||
import h5py
|
||||
self._f = h5py.File(filename, 'r')
|
||||
|
||||
# Ensure filetype and revision are correct
|
||||
if 'filetype' not in self._f or self._f[
|
||||
'filetype'].value.decode() != 'statepoint':
|
||||
raise IOError('{} is not a statepoint file.'.format(filename))
|
||||
if self._f['revision'].value != 14:
|
||||
raise IOError('Statepoint file has a file revision of {} '
|
||||
'which is not consistent with the revision this '
|
||||
'version of OpenMC expects ({}).'.format(
|
||||
self._f['revision'].value, 14))
|
||||
|
||||
# Set flags for what data has been read
|
||||
self._results = False
|
||||
self._source = False
|
||||
self._meshes_read = False
|
||||
self._tallies_read = False
|
||||
self._with_summary = False
|
||||
|
||||
# Read all metadata
|
||||
self._read_metadata()
|
||||
|
||||
# Read information about tally meshes
|
||||
self._read_meshes()
|
||||
|
||||
# Read tally metadata
|
||||
self._read_tallies()
|
||||
self._global_tallies = None
|
||||
|
||||
def close(self):
|
||||
self._f.close()
|
||||
|
||||
@property
|
||||
def k_combined(self):
|
||||
return self._k_combined
|
||||
def cmfd_on(self):
|
||||
return self._f['cmfd_on'].value > 0
|
||||
|
||||
@property
|
||||
def n_particles(self):
|
||||
return self._n_particles
|
||||
def cmfd_balance(self):
|
||||
if self.cmfd_on:
|
||||
return self._f['cmfd/cmfd_balance'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def n_batches(self):
|
||||
return self._n_batches
|
||||
def cmfd_dominance(self):
|
||||
if self.cmfd_on:
|
||||
return self._f['cmfd/cmfd_dominance'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def cmfd_entropy(self):
|
||||
if self.cmfd_on:
|
||||
return self._f['cmfd/cmfd_entropy'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def cmfd_indices(self):
|
||||
if self.cmfd_on:
|
||||
return self._f['cmfd/indices'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def cmfd_src(self):
|
||||
if self.cmfd_on:
|
||||
data = self._f['cmfd/cmfd_src'].value
|
||||
return np.reshape(data, tuple(self.cmfd_indices), order='F')
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def cmfd_srccmp(self):
|
||||
if self.cmfd_on:
|
||||
return self._f['cmfd/cmfd_srccmp'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def current_batch(self):
|
||||
return self._current_batch
|
||||
return self._f['current_batch'].value
|
||||
|
||||
@property
|
||||
def results(self):
|
||||
return self._results
|
||||
def date_and_time(self):
|
||||
return self._f['date_and_time'].value.decode()
|
||||
|
||||
@property
|
||||
def entropy(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['entropy'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def gen_per_batch(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['gen_per_batch'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def global_tallies(self):
|
||||
if self._global_tallies is None:
|
||||
data = self._f['global_tallies'].value
|
||||
gt = np.zeros_like(data, dtype=[
|
||||
('name', 'a14'), ('sum', 'f8'), ('sum_sq', 'f8'),
|
||||
('mean', 'f8'), ('std_dev', 'f8')])
|
||||
gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength',
|
||||
'leakage']
|
||||
gt['sum'] = data['sum']
|
||||
gt['sum_sq'] = data['sum_sq']
|
||||
|
||||
# Calculate mean and sample standard deviation of mean
|
||||
n = self.n_realizations
|
||||
gt['mean'] = gt['sum']/n
|
||||
gt['std_dev'] = np.sqrt((gt['sum_sq']/n - gt['mean']**2)/(n - 1))
|
||||
|
||||
self._global_tallies = gt
|
||||
|
||||
return self._global_tallies
|
||||
|
||||
@property
|
||||
def k_cmfd(self):
|
||||
if self.cmfd_on:
|
||||
return self._f['cmfd/k_cmfd'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def k_generation(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['k_generation'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def k_combined(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['k_combined'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def k_col_abs(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['k_col_abs'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def k_col_tra(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['k_col_tra'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def k_abs_tra(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['k_abs_tra'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def meshes(self):
|
||||
if not self._meshes_read:
|
||||
# Initialize dictionaries for the Meshes
|
||||
# Keys - Mesh IDs
|
||||
# Values - Mesh objects
|
||||
self._meshes = {}
|
||||
|
||||
# Read the number of Meshes
|
||||
n_meshes = self._f['tallies/meshes/n_meshes'].value
|
||||
|
||||
# Read a list of the IDs for each Mesh
|
||||
if n_meshes > 0:
|
||||
# User-defined Mesh IDs
|
||||
mesh_keys = self._f['tallies/meshes/keys'].value
|
||||
else:
|
||||
mesh_keys = []
|
||||
|
||||
# Build dictionary of Meshes
|
||||
base = 'tallies/meshes/mesh '
|
||||
|
||||
# Iterate over all Meshes
|
||||
for mesh_key in mesh_keys:
|
||||
# Read the mesh type
|
||||
mesh_type = self._f['{0}{1}/type'.format(base, mesh_key)].value.decode()
|
||||
|
||||
# Read the mesh dimensions, lower-left coordinates,
|
||||
# upper-right coordinates, and width of each mesh cell
|
||||
dimension = self._f['{0}{1}/dimension'.format(base, mesh_key)].value
|
||||
lower_left = self._f['{0}{1}/lower_left'.format(base, mesh_key)].value
|
||||
upper_right = self._f['{0}{1}/upper_right'.format(base, mesh_key)].value
|
||||
width = self._f['{0}{1}/width'.format(base, mesh_key)].value
|
||||
|
||||
# Create the Mesh and assign properties to it
|
||||
mesh = openmc.Mesh(mesh_key)
|
||||
mesh.dimension = dimension
|
||||
mesh.width = width
|
||||
mesh.lower_left = lower_left
|
||||
mesh.upper_right = upper_right
|
||||
mesh.type = mesh_type
|
||||
|
||||
# Add mesh to the global dictionary of all Meshes
|
||||
self._meshes[mesh_key] = mesh
|
||||
|
||||
self._meshes_read = True
|
||||
|
||||
return self._meshes
|
||||
|
||||
@property
|
||||
def n_batches(self):
|
||||
return self._f['n_batches'].value
|
||||
|
||||
@property
|
||||
def n_inactive(self):
|
||||
if self.run_mode == 'k-eigenvalue':
|
||||
return self._f['n_inactive'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def n_particles(self):
|
||||
return self._f['n_particles'].value
|
||||
|
||||
@property
|
||||
def n_realizations(self):
|
||||
return self._f['n_realizations'].value
|
||||
|
||||
@property
|
||||
def path(self):
|
||||
return self._f['path'].value.decode()
|
||||
|
||||
@property
|
||||
def run_mode(self):
|
||||
return self._f['run_mode'].value.decode()
|
||||
|
||||
@property
|
||||
def seed(self):
|
||||
return self._f['seed'].value
|
||||
|
||||
@property
|
||||
def source(self):
|
||||
return self._source
|
||||
if self.source_present:
|
||||
return self._f['source_bank'].value
|
||||
else:
|
||||
return None
|
||||
|
||||
@property
|
||||
def with_summary(self):
|
||||
return self._with_summary
|
||||
def source_present(self):
|
||||
return self._f['source_present'].value > 0
|
||||
|
||||
@property
|
||||
def tallies(self):
|
||||
if not self._tallies_read:
|
||||
# Initialize dictionary for tallies
|
||||
self._tallies = {}
|
||||
|
||||
# Read the number of tallies
|
||||
n_tallies = self._f['tallies/n_tallies'].value
|
||||
|
||||
# Read a list of the IDs for each Tally
|
||||
if n_tallies > 0:
|
||||
# OpenMC Tally IDs (redefined internally from user definitions)
|
||||
tally_keys = self._f['tallies/keys'].value
|
||||
else:
|
||||
tally_keys = []
|
||||
|
||||
base = 'tallies/tally '
|
||||
|
||||
# Iterate over all Tallies
|
||||
for tally_key in tally_keys:
|
||||
|
||||
# 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_id=tally_key)
|
||||
tally._statepoint = self
|
||||
tally.estimator = self._f['{0}{1}/estimator'.format(
|
||||
base, tally_key)].value.decode()
|
||||
tally.num_realizations = n_realizations
|
||||
|
||||
# Read the number of Filters
|
||||
n_filters = self._f['{0}{1}/n_filters'.format(base, tally_key)].value
|
||||
|
||||
subbase = '{0}{1}/filter '.format(base, tally_key)
|
||||
|
||||
# Initialize all Filters
|
||||
for j in range(1, n_filters+1):
|
||||
|
||||
# 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
|
||||
|
||||
n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value
|
||||
|
||||
# Read the bin values
|
||||
bins = self._f['{0}{1}/bins'.format(subbase, j)].value
|
||||
|
||||
# Create Filter object
|
||||
filter = openmc.Filter(filter_type, bins)
|
||||
filter.offset = offset
|
||||
filter.num_bins = n_bins
|
||||
|
||||
if filter_type == 'mesh':
|
||||
mesh_ids = self._f['tallies/meshes/ids'].value
|
||||
mesh_keys = self._f['tallies/meshes/keys'].value
|
||||
|
||||
key = mesh_keys[mesh_ids == bins][0]
|
||||
filter.mesh = self.meshes[key]
|
||||
|
||||
# Add Filter to the Tally
|
||||
tally.add_filter(filter)
|
||||
|
||||
# Read Nuclide bins
|
||||
nuclide_names = self._f['{0}{1}/nuclides'.format(base, tally_key)].value
|
||||
|
||||
# Add all Nuclides to the Tally
|
||||
for name in nuclide_names:
|
||||
nuclide = openmc.Nuclide(name.decode().strip())
|
||||
tally.add_nuclide(nuclide)
|
||||
|
||||
# Read score bins
|
||||
n_score_bins = self._f['{0}{1}/n_score_bins'.format(base, tally_key)].value
|
||||
|
||||
tally.num_score_bins = n_score_bins
|
||||
|
||||
scores = self._f['{0}{1}/score_bins'.format(
|
||||
base, tally_key)].value
|
||||
n_user_scores = self._f['{0}{1}/n_user_score_bins'
|
||||
.format(base, tally_key)].value
|
||||
|
||||
# Compute and set the filter strides
|
||||
for i in range(n_filters):
|
||||
filter = tally.filters[i]
|
||||
filter.stride = n_score_bins * len(nuclide_names)
|
||||
|
||||
for j in range(i+1, n_filters):
|
||||
filter.stride *= tally.filters[j].num_bins
|
||||
|
||||
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
|
||||
moments = self._f['{0}{1}/moment_orders'.format(
|
||||
base, tally_key)].value
|
||||
|
||||
# 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())
|
||||
elif '-pn' in score:
|
||||
score = score.replace('-pn', '-' + moments[j].decode())
|
||||
elif '-yn' in score:
|
||||
score = score.replace('-yn', '-' + moments[j].decode())
|
||||
|
||||
tally.add_score(score)
|
||||
|
||||
# Add Tally to the global dictionary of all Tallies
|
||||
self._tallies[tally_key] = tally
|
||||
|
||||
self._tallies_read = True
|
||||
|
||||
return self._tallies
|
||||
|
||||
@property
|
||||
def tallies_present(self):
|
||||
return self._tallies_present
|
||||
return self._f['tallies/tallies_present'].value
|
||||
|
||||
@property
|
||||
def global_tallies(self):
|
||||
return self._global_tallies
|
||||
def version(self):
|
||||
return (self._f['version_major'].value,
|
||||
self._f['version_minor'].value,
|
||||
self._f['version_release'].value)
|
||||
|
||||
@property
|
||||
def n_realizations(self):
|
||||
return self._n_realizations
|
||||
|
||||
def _read_metadata(self):
|
||||
# Read filetype
|
||||
self._filetype = self._f['filetype'].value
|
||||
|
||||
# Read statepoint revision
|
||||
self._revision = self._f['revision'].value
|
||||
if self._revision != 13:
|
||||
raise Exception('Statepoint Revision is not consistent.')
|
||||
|
||||
# Read OpenMC version
|
||||
self._version = [self._f['version_major'].value,
|
||||
self._f['version_minor'].value,
|
||||
self._f['version_release'].value]
|
||||
|
||||
# Read date and time
|
||||
self._date_and_time = self._f['date_and_time'].value[0]
|
||||
|
||||
# Read path
|
||||
self._path = self._f['path'].value[0].strip()
|
||||
|
||||
# Read random number seed
|
||||
self._seed = self._f['seed'].value
|
||||
|
||||
# Read run information
|
||||
self._run_mode = self._f['run_mode'].value
|
||||
self._n_particles = self._f['n_particles'].value
|
||||
self._n_batches = self._f['n_batches'].value
|
||||
|
||||
# Read current batch
|
||||
self._current_batch = self._f['current_batch'].value
|
||||
|
||||
# Read whether or not the source site distribution is present
|
||||
self._source_present = self._f['source_present'].value
|
||||
|
||||
# Read criticality information
|
||||
if self._run_mode == 2:
|
||||
self._read_criticality()
|
||||
|
||||
def _read_criticality(self):
|
||||
# Read criticality information
|
||||
if self._run_mode == 2:
|
||||
|
||||
self._n_inactive = self._f['n_inactive'].value
|
||||
self._gen_per_batch = self._f['gen_per_batch'].value
|
||||
self._k_generation = self._f['k_generation'].value
|
||||
self._entropy = self._f['entropy'].value
|
||||
|
||||
self._k_col_abs = self._f['k_col_abs'].value
|
||||
self._k_col_tra = self._f['k_col_tra'].value
|
||||
self._k_abs_tra = self._f['k_abs_tra'].value
|
||||
self._k_combined = self._f['k_combined'].value
|
||||
|
||||
# Read CMFD information (if used)
|
||||
self._read_cmfd()
|
||||
|
||||
def _read_cmfd(self):
|
||||
base = 'cmfd'
|
||||
|
||||
# Read CMFD information
|
||||
self._cmfd_on = self._f['cmfd_on'].value
|
||||
|
||||
if self._cmfd_on == 1:
|
||||
self._cmfd_indices = self._f['{0}/indices'.format(base)].value
|
||||
self._k_cmfd = self._f['{0}/k_cmfd'.format(base)].value
|
||||
self._cmfd_src = self._f['{0}/cmfd_src'.format(base)].value
|
||||
self._cmfd_src = np.reshape(self._cmfd_src, tuple(self._cmfd_indices),
|
||||
order='F')
|
||||
self._cmfd_entropy = self._f['{0}/cmfd_entropy'.format(base)].value
|
||||
self._cmfd_balance = self._f['{0}/cmfd_balance'.format(base)].value
|
||||
self._cmfd_dominance = self._f['{0}/cmfd_dominance'.format(base)].value
|
||||
self._cmfd_srccmp = self._f['{0}/cmfd_srccmp'.format(base)].value
|
||||
|
||||
def _read_meshes(self):
|
||||
# Initialize dictionaries for the Meshes
|
||||
# Keys - Mesh IDs
|
||||
# Values - Mesh objects
|
||||
self._meshes = {}
|
||||
|
||||
# Read the number of Meshes
|
||||
self._n_meshes = self._f['tallies/meshes/n_meshes'].value
|
||||
|
||||
# Read a list of the IDs for each Mesh
|
||||
if self._n_meshes > 0:
|
||||
|
||||
# OpenMC Mesh IDs (redefined internally from user definitions)
|
||||
self._mesh_ids = self._f['tallies/meshes/ids'].value
|
||||
|
||||
# User-defined Mesh IDs
|
||||
self._mesh_keys = self._f['tallies/meshes/keys'].value
|
||||
|
||||
else:
|
||||
self._mesh_keys = []
|
||||
self._mesh_ids = []
|
||||
|
||||
# Build dictionary of Meshes
|
||||
base = 'tallies/meshes/mesh '
|
||||
|
||||
# Iterate over all Meshes
|
||||
for mesh_key in self._mesh_keys:
|
||||
|
||||
# Read the user-specified Mesh ID and type
|
||||
mesh_id = self._f['{0}{1}/id'.format(base, mesh_key)].value
|
||||
mesh_type = self._f['{0}{1}/type'.format(base, mesh_key)].value
|
||||
|
||||
# Get the Mesh dimension
|
||||
n_dimension = self._f['{0}{1}/n_dimension'.format(base, mesh_key)].value
|
||||
|
||||
# Read the mesh dimensions, lower-left coordinates,
|
||||
# upper-right coordinates, and width of each mesh cell
|
||||
dimension = self._f['{0}{1}/dimension'.format(base, mesh_key)].value
|
||||
lower_left = self._f['{0}{1}/lower_left'.format(base, mesh_key)].value
|
||||
upper_right = self._f['{0}{1}/upper_right'.format(base, mesh_key)].value
|
||||
width = self._f['{0}{1}/width'.format(base, mesh_key)].value
|
||||
|
||||
# Create the Mesh and assign properties to it
|
||||
mesh = openmc.Mesh(mesh_id)
|
||||
|
||||
mesh.dimension = dimension
|
||||
mesh.width = width
|
||||
mesh.lower_left = lower_left
|
||||
mesh.upper_right = upper_right
|
||||
|
||||
#FIXME: Set the mesh type to 'rectangular' by default
|
||||
mesh.type = 'rectangular'
|
||||
|
||||
# Add mesh to the global dictionary of all Meshes
|
||||
self._meshes[mesh_id] = mesh
|
||||
|
||||
def _read_tallies(self):
|
||||
# Initialize dictionaries for the Tallies
|
||||
# Keys - Tally IDs
|
||||
# Values - Tally objects
|
||||
self._tallies = {}
|
||||
|
||||
# Read the number of tallies
|
||||
self._n_tallies = self._f['/tallies/n_tallies'].value
|
||||
|
||||
# Read a list of the IDs for each Tally
|
||||
if self._n_tallies > 0:
|
||||
|
||||
# OpenMC Tally IDs (redefined internally from user definitions)
|
||||
self._tally_ids = self._f['tallies/ids'].value
|
||||
|
||||
# User-defined Tally IDs
|
||||
self._tally_keys = self._f['tallies/keys'].value
|
||||
|
||||
else:
|
||||
self._tally_keys = []
|
||||
self._tally_ids = []
|
||||
|
||||
base = 'tallies/tally '
|
||||
|
||||
# Iterate over all Tallies
|
||||
for tally_key in self._tally_keys:
|
||||
|
||||
# Read integer Tally estimator type code (analog or tracklength)
|
||||
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.estimator = ESTIMATOR_TYPES[estimator_type]
|
||||
tally.num_realizations = n_realizations
|
||||
|
||||
# Read the number of Filters
|
||||
n_filters = self._f['{0}{1}/n_filters'.format(base, tally_key)].value
|
||||
|
||||
subbase = '{0}{1}/filter '.format(base, tally_key)
|
||||
|
||||
# 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 offset
|
||||
offset = self._f['{0}{1}/offset'.format(subbase, j)].value
|
||||
|
||||
n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value
|
||||
|
||||
if n_bins <= 0:
|
||||
msg = 'Unable to create Filter "{0}" for Tally ID="{1}" ' \
|
||||
'since no bins were specified'.format(j, tally_key)
|
||||
raise ValueError(msg)
|
||||
|
||||
# Read the bin values
|
||||
if FILTER_TYPES[filter_type] in ['energy', 'energyout']:
|
||||
bins = self._f['{0}{1}/bins'.format(subbase, j)].value
|
||||
|
||||
elif FILTER_TYPES[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.offset = offset
|
||||
filter.num_bins = n_bins
|
||||
|
||||
if FILTER_TYPES[filter_type] == 'mesh':
|
||||
key = self._mesh_keys[self._mesh_ids == bins][0]
|
||||
filter.mesh = self._meshes[key]
|
||||
|
||||
# Add Filter to the Tally
|
||||
tally.add_filter(filter)
|
||||
|
||||
# Read Nuclide bins
|
||||
n_nuclides = self._f['{0}{1}/n_nuclides'.format(base, tally_key)].value
|
||||
|
||||
nuclide_zaids = self._f['{0}{1}/nuclides'.format(base, tally_key)].value
|
||||
|
||||
# Add all Nuclides to the Tally
|
||||
for nuclide_zaid in nuclide_zaids:
|
||||
tally.add_nuclide(nuclide_zaid)
|
||||
|
||||
# Read score bins
|
||||
n_score_bins = self._f['{0}{1}/n_score_bins'.format(base, tally_key)].value
|
||||
|
||||
tally.num_score_bins = n_score_bins
|
||||
|
||||
score_bins = self._f['{0}{1}/score_bins'.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
|
||||
|
||||
# Compute and set the filter strides
|
||||
for i in range(n_filters):
|
||||
filter = tally.filters[i]
|
||||
filter.stride = n_score_bins * n_nuclides
|
||||
|
||||
for j in range(i+1, n_filters):
|
||||
filter.stride *= tally.filters[j].num_bins
|
||||
|
||||
# Read scattering moment order strings (e.g., P3, Y-1,2, etc.)
|
||||
moments = []
|
||||
subbase = '{0}{1}/moments/'.format(base, tally_key)
|
||||
|
||||
# Extract the moment order string for each score
|
||||
for k in range(len(scores)):
|
||||
moment = str(self._f['{0}order{1}'.format(
|
||||
subbase, k+1)].value[0])
|
||||
moment = moment.lstrip('[\'')
|
||||
moment = moment.rstrip('\']')
|
||||
|
||||
# Remove extra whitespace
|
||||
moment.replace(" ", "")
|
||||
moments.append(moment)
|
||||
|
||||
# Add the scores to the Tally
|
||||
for j, score in enumerate(scores):
|
||||
# If this is a scattering moment, insert the scattering order
|
||||
if '-n' in score:
|
||||
score = score.replace('-n', '-' + str(moments[j]))
|
||||
elif '-pn' in score:
|
||||
score = score.replace('-pn', '-' + str(moments[j]))
|
||||
elif '-yn' in score:
|
||||
score = score.replace('-yn', '-' + str(moments[j]))
|
||||
|
||||
tally.add_score(score)
|
||||
|
||||
# Add Tally to the global dictionary of all Tallies
|
||||
self.tallies[tally_key] = tally
|
||||
|
||||
def read_results(self):
|
||||
"""Read tally results and store them in the ``tallies`` attribute. No results
|
||||
are read when the statepoint is instantiated.
|
||||
|
||||
"""
|
||||
|
||||
# Number of realizations for global Tallies
|
||||
self._n_realizations = self._f['n_realizations'].value
|
||||
|
||||
# Read global Tallies
|
||||
n_global_tallies = self._f['n_global_tallies'].value
|
||||
|
||||
data = self._f['global_tallies'].value
|
||||
self._global_tallies = np.column_stack((data['sum'], data['sum_sq']))
|
||||
|
||||
# Flag indicating if Tallies are present
|
||||
self._tallies_present = self._f['tallies/tallies_present'].value
|
||||
|
||||
base = 'tallies/tally '
|
||||
|
||||
# Read Tally results
|
||||
if self._tallies_present:
|
||||
|
||||
# Iterate over and extract the results for all Tallies
|
||||
for tally_key in self._tally_keys:
|
||||
|
||||
# Get this Tally
|
||||
tally = self._tallies[tally_key]
|
||||
|
||||
# Compute the total number of bins for this Tally
|
||||
num_tot_bins = tally.num_bins
|
||||
|
||||
# Extract Tally data from the file
|
||||
data = self._f['{0}{1}/results'.format(base, tally_key)].value
|
||||
sum = data['sum']
|
||||
sum_sq = data['sum_sq']
|
||||
|
||||
# Define a routine to convert 0 to 1
|
||||
def nonzero(val):
|
||||
return 1 if not val else val
|
||||
|
||||
# Reshape the results arrays
|
||||
new_shape = (nonzero(tally.num_filter_bins),
|
||||
nonzero(tally.num_nuclides),
|
||||
nonzero(tally.num_score_bins))
|
||||
sum = np.reshape(sum, new_shape)
|
||||
sum_sq = np.reshape(sum_sq, new_shape)
|
||||
|
||||
# Set the data for this Tally
|
||||
tally.sum = sum
|
||||
tally.sum_sq = sum_sq
|
||||
|
||||
# Indicate that Tally results have been read
|
||||
self._results = True
|
||||
|
||||
def read_source(self):
|
||||
"""Read and store source sites from the statepoint file. By default, source
|
||||
sites are not loaded upon initialization.
|
||||
|
||||
"""
|
||||
|
||||
# Check whether Tally results have been read
|
||||
if not self._results:
|
||||
self.read_results()
|
||||
|
||||
# Check if source bank is in statepoint
|
||||
if not self._source_present:
|
||||
print('Unable to read source since it is not in statepoint file')
|
||||
return
|
||||
|
||||
# Initialize a NumPy array for the source sites
|
||||
self._source = np.empty(self._n_particles, dtype=SourceSite)
|
||||
|
||||
# For HDF5 state points, copy entire bank
|
||||
source_sites = self._f['source_bank'].value
|
||||
|
||||
# Initialize SourceSite object for each particle
|
||||
for i in range(self._n_particles):
|
||||
# Initialize new source site
|
||||
site = SourceSite()
|
||||
|
||||
# Read position, angle, and energy
|
||||
site._weight, site._xyz, site._uvw, site._E = source_sites[i]
|
||||
|
||||
# Store the source site in the NumPy array
|
||||
self._source[i] = site
|
||||
|
||||
def compute_ci(self, confidence=0.95):
|
||||
"""Computes confidence intervals for each Tally bin.
|
||||
|
||||
This method is equivalent to calling compute_stdev(...) when the
|
||||
confidence is known as opposed to its corresponding t value.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
confidence : float, optional
|
||||
Confidence level. Defaults to 0.95.
|
||||
|
||||
"""
|
||||
|
||||
# Determine significance level and percentile for two-sided CI
|
||||
alpha = 1 - confidence
|
||||
percentile = 1 - alpha/2
|
||||
|
||||
# Calculate t-value
|
||||
t_value = scipy.stats.t.ppf(percentile, self._n_realizations - 1)
|
||||
self.compute_stdev(t_value)
|
||||
|
||||
def compute_stdev(self, t_value=1.0):
|
||||
"""Computes the sample mean and the standard deviation of the mean
|
||||
for each Tally bin.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
t_value : float, optional
|
||||
Student's t-value applied to the uncertainty. Defaults to 1.0,
|
||||
meaning the reported value is the sample standard deviation.
|
||||
|
||||
"""
|
||||
|
||||
# Determine number of realizations
|
||||
n = self._n_realizations
|
||||
|
||||
# Calculate the standard deviation for each global tally
|
||||
for i in range(len(self._global_tallies)):
|
||||
|
||||
# Get sum and sum of squares
|
||||
s, s2 = self._global_tallies[i]
|
||||
|
||||
# Calculate sample mean and replace value
|
||||
s /= n
|
||||
self._global_tallies[i, 0] = s
|
||||
|
||||
# Calculate standard deviation
|
||||
if s != 0.0:
|
||||
self._global_tallies[i, 1] = t_value * np.sqrt((s2 / n - s**2) / (n-1))
|
||||
|
||||
# Calculate sample mean and standard deviation for user-defined Tallies
|
||||
for tally_id, tally in self.tallies.items():
|
||||
tally.compute_std_dev(t_value)
|
||||
def with_summary(self):
|
||||
return self._with_summary
|
||||
|
||||
def get_tally(self, scores=[], filters=[], nuclides=[],
|
||||
name=None, id=None, estimator=None):
|
||||
|
|
@ -706,15 +604,6 @@ class StatePoint(object):
|
|||
tally.name = summary.tallies[tally_id].name
|
||||
tally.with_summary = True
|
||||
|
||||
nuclide_zaids = copy.deepcopy(tally.nuclides)
|
||||
|
||||
for nuclide_zaid in nuclide_zaids:
|
||||
tally.remove_nuclide(nuclide_zaid)
|
||||
if nuclide_zaid == -1:
|
||||
tally.add_nuclide(openmc.Nuclide('total'))
|
||||
else:
|
||||
tally.add_nuclide(summary.nuclides[nuclide_zaid])
|
||||
|
||||
for filter in tally.filters:
|
||||
if filter.type == 'surface':
|
||||
surface_ids = []
|
||||
|
|
|
|||
|
|
@ -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'][0]
|
||||
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
|
||||
|
||||
|
|
@ -97,23 +67,17 @@ class Summary(object):
|
|||
|
||||
material_id = int(key.lstrip('material '))
|
||||
index = self._f['materials'][key]['index'].value
|
||||
name = self._f['materials'][key]['name'][0]
|
||||
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
|
||||
|
||||
sab_names = []
|
||||
sab_xs = []
|
||||
|
||||
# Read the names of the S(a,b) tables for this Material
|
||||
for i in range(1, n_sab+1):
|
||||
sab_table = \
|
||||
self._f['materials'][key]['sab_tables'][str(i)].value
|
||||
|
||||
# Read the cross-section identifiers for each S(a,b) table
|
||||
sab_names.append(sab_table.split('.')[0])
|
||||
sab_xs.append(sab_table.split('.')[1])
|
||||
# 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
|
||||
for sab_table in sab_tables:
|
||||
name, xs = sab_table.decode().split('.')
|
||||
material.add_s_alpha_beta(name, xs)
|
||||
|
||||
# Create the Material
|
||||
material = openmc.Material(material_id=material_id, name=name)
|
||||
|
|
@ -121,21 +85,17 @@ 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 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)
|
||||
if 'nat' in name:
|
||||
material.add_element(openmc.Element(name=name, xs=xs),
|
||||
percent=density, percent_type='ao')
|
||||
else:
|
||||
material.add_nuclide(openmc.Nuclide(name=name, xs=xs),
|
||||
percent=density, percent_type='ao')
|
||||
|
||||
# Add the Material to the global dictionary of all Materials
|
||||
self.materials[index] = material
|
||||
|
|
@ -154,67 +114,67 @@ class Summary(object):
|
|||
|
||||
surface_id = int(key.lstrip('surface '))
|
||||
index = self._f['geometry/surfaces'][key]['index'].value
|
||||
name = self._f['geometry/surfaces'][key]['name'][0]
|
||||
surf_type = self._f['geometry/surfaces'][key]['type'][...]
|
||||
bc = self._f['geometry/surfaces'][key]['boundary_condition'][...][0]
|
||||
name = self._f['geometry/surfaces'][key]['name'].value.decode()
|
||||
surf_type = self._f['geometry/surfaces'][key]['type'].value.decode()
|
||||
bc = self._f['geometry/surfaces'][key]['boundary_condition'].value.decode()
|
||||
coeffs = self._f['geometry/surfaces'][key]['coefficients'][...]
|
||||
|
||||
# Create the Surface based on its type
|
||||
if surf_type == 'X Plane':
|
||||
if surf_type == 'x-plane':
|
||||
x0 = coeffs[0]
|
||||
surface = openmc.XPlane(surface_id, bc, x0, name)
|
||||
|
||||
elif surf_type == 'Y Plane':
|
||||
elif surf_type == 'y-plane':
|
||||
y0 = coeffs[0]
|
||||
surface = openmc.YPlane(surface_id, bc, y0, name)
|
||||
|
||||
elif surf_type == 'Z Plane':
|
||||
elif surf_type == 'z-plane':
|
||||
z0 = coeffs[0]
|
||||
surface = openmc.ZPlane(surface_id, bc, z0, name)
|
||||
|
||||
elif surf_type == 'Plane':
|
||||
elif surf_type == 'plane':
|
||||
A = coeffs[0]
|
||||
B = coeffs[1]
|
||||
C = coeffs[2]
|
||||
D = coeffs[3]
|
||||
surface = openmc.Plane(surface_id, bc, A, B, C, D, name)
|
||||
|
||||
elif surf_type == 'X Cylinder':
|
||||
elif surf_type == 'x-cylinder':
|
||||
y0 = coeffs[0]
|
||||
z0 = coeffs[1]
|
||||
R = coeffs[2]
|
||||
surface = openmc.XCylinder(surface_id, bc, y0, z0, R, name)
|
||||
|
||||
elif surf_type == 'Y Cylinder':
|
||||
elif surf_type == 'y-cylinder':
|
||||
x0 = coeffs[0]
|
||||
z0 = coeffs[1]
|
||||
R = coeffs[2]
|
||||
surface = openmc.YCylinder(surface_id, bc, x0, z0, R, name)
|
||||
|
||||
elif surf_type == 'Z Cylinder':
|
||||
elif surf_type == 'z-cylinder':
|
||||
x0 = coeffs[0]
|
||||
y0 = coeffs[1]
|
||||
R = coeffs[2]
|
||||
surface = openmc.ZCylinder(surface_id, bc, x0, y0, R, name)
|
||||
|
||||
elif surf_type == 'Sphere':
|
||||
elif surf_type == 'sphere':
|
||||
x0 = coeffs[0]
|
||||
y0 = coeffs[1]
|
||||
z0 = coeffs[2]
|
||||
R = coeffs[3]
|
||||
surface = openmc.Sphere(surface_id, bc, x0, y0, z0, R, name)
|
||||
|
||||
elif surf_type in ['X Cone', 'Y Cone', 'Z Cone']:
|
||||
elif surf_type in ['x-cone', 'y-cone', 'z-cone']:
|
||||
x0 = coeffs[0]
|
||||
y0 = coeffs[1]
|
||||
z0 = coeffs[2]
|
||||
R2 = coeffs[3]
|
||||
|
||||
if surf_type == 'X Cone':
|
||||
if surf_type == 'x-cone':
|
||||
surface = openmc.XCone(surface_id, bc, x0, y0, z0, R2, name)
|
||||
if surf_type == 'Y Cone':
|
||||
if surf_type == 'y-cone':
|
||||
surface = openmc.YCone(surface_id, bc, x0, y0, z0, R2, name)
|
||||
if surf_type == 'Z Cone':
|
||||
if surf_type == 'z-cone':
|
||||
surface = openmc.ZCone(surface_id, bc, x0, y0, z0, R2, name)
|
||||
|
||||
# Add Surface to global dictionary of all Surfaces
|
||||
|
|
@ -242,8 +202,8 @@ class Summary(object):
|
|||
|
||||
cell_id = int(key.lstrip('cell '))
|
||||
index = self._f['geometry/cells'][key]['index'].value
|
||||
name = self._f['geometry/cells'][key]['name'][0]
|
||||
fill_type = self._f['geometry/cells'][key]['fill_type'][...][0]
|
||||
name = self._f['geometry/cells'][key]['name'].value.decode()
|
||||
fill_type = self._f['geometry/cells'][key]['fill_type'].value.decode()
|
||||
|
||||
if fill_type == 'normal':
|
||||
fill = self._f['geometry/cells'][key]['material'].value
|
||||
|
|
@ -261,21 +221,17 @@ class Summary(object):
|
|||
cell = openmc.Cell(cell_id=cell_id, name=name)
|
||||
|
||||
if fill_type == 'universe':
|
||||
maps = self._f['geometry/cells'][key]['maps'].value
|
||||
|
||||
if maps > 0:
|
||||
if 'offset' in self._f['geometry/cells'][key]:
|
||||
offset = self._f['geometry/cells'][key]['offset'][...]
|
||||
cell.offsets = offset
|
||||
|
||||
translated = self._f['geometry/cells'][key]['translated'].value
|
||||
if translated:
|
||||
if 'translation' in self._f['geometry/cells'][key]:
|
||||
translation = \
|
||||
self._f['geometry/cells'][key]['translation'][...]
|
||||
translation = np.asarray(translation, dtype=np.float64)
|
||||
cell.translation = translation
|
||||
|
||||
rotated = self._f['geometry/cells'][key]['rotated'].value
|
||||
if rotated:
|
||||
if 'rotation' in self._f['geometry/cells'][key]:
|
||||
rotation = \
|
||||
self._f['geometry/cells'][key]['rotation'][...]
|
||||
rotation = np.asarray(rotation, dtype=np.int)
|
||||
|
|
@ -288,8 +244,8 @@ class Summary(object):
|
|||
for surface_halfspace in surfaces:
|
||||
|
||||
halfspace = np.sign(surface_halfspace)
|
||||
surface_id = np.abs(surface_halfspace)
|
||||
surface = self.surfaces[surface_id]
|
||||
surface_id = abs(surface_halfspace)
|
||||
surface = self.get_surface_by_id(surface_id)
|
||||
cell.add_surface(surface, halfspace)
|
||||
|
||||
# Add the Cell to the global dictionary of all Cells
|
||||
|
|
@ -336,13 +292,13 @@ class Summary(object):
|
|||
|
||||
lattice_id = int(key.lstrip('lattice '))
|
||||
index = self._f['geometry/lattices'][key]['index'].value
|
||||
name = self._f['geometry/lattices'][key]['name'][...][0]
|
||||
lattice_type = self._f['geometry/lattices'][key]['type'][...][0]
|
||||
maps = self._f['geometry/lattices'][key]['maps'].value
|
||||
offset_size = self._f['geometry/lattices'][key]['offset_size'].value
|
||||
name = self._f['geometry/lattices'][key]['name'].value.decode()
|
||||
lattice_type = self._f['geometry/lattices'][key]['type'].value.decode()
|
||||
|
||||
if offset_size > 0:
|
||||
if 'offsets' in self._f['geometry/lattices'][key]:
|
||||
offsets = self._f['geometry/lattices'][key]['offsets'][...]
|
||||
else:
|
||||
offsets = None
|
||||
|
||||
if lattice_type == 'rectangular':
|
||||
dimension = self._f['geometry/lattices'][key]['dimension'][...]
|
||||
|
|
@ -383,7 +339,7 @@ class Summary(object):
|
|||
universes = universes[:, ::-1, :]
|
||||
lattice.universes = universes
|
||||
|
||||
if offset_size > 0:
|
||||
if offsets is not None:
|
||||
offsets = np.swapaxes(offsets, 0, 1)
|
||||
offsets = np.swapaxes(offsets, 1, 2)
|
||||
lattice.offsets = offsets
|
||||
|
|
@ -477,7 +433,7 @@ class Summary(object):
|
|||
# Lattice is 2D; extract the only axial level
|
||||
lattice.universes = universes[0]
|
||||
|
||||
if offset_size > 0:
|
||||
if offsets is not None:
|
||||
lattice.offsets = offsets
|
||||
|
||||
# Add the Lattice to the global dictionary of all Lattices
|
||||
|
|
@ -532,26 +488,19 @@ class Summary(object):
|
|||
|
||||
# Iterate over all Tallies
|
||||
for tally_key in tally_keys:
|
||||
|
||||
tally_id = int(tally_key.strip('tally '))
|
||||
subbase = '{0}{1}'.format(base, tally_id)
|
||||
|
||||
# Read Tally name metadata
|
||||
name_size = self._f['{0}/name_size'.format(subbase)][...]
|
||||
if (name_size > 0):
|
||||
tally_name = self._f['{0}/name'.format(subbase)][...][0]
|
||||
tally_name = tally_name.lstrip('[\'')
|
||||
tally_name = tally_name.rstrip('\']')
|
||||
else:
|
||||
tally_name = ''
|
||||
tally_name = self._f['{0}/name'.format(subbase)].value.decode()
|
||||
|
||||
# Create Tally object and assign basic properties
|
||||
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}/score_bins'.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
|
||||
|
||||
|
|
@ -560,12 +509,10 @@ class Summary(object):
|
|||
|
||||
# Initialize all Filters
|
||||
for j in range(1, num_filters+1):
|
||||
|
||||
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
|
||||
|
|
@ -597,29 +544,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
|
||||
|
||||
|
|
|
|||
|
|
@ -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):
|
||||
|
|
@ -134,7 +135,7 @@ class Surface(object):
|
|||
|
||||
element.set("type", self._type)
|
||||
element.set("boundary", self._boundary_type)
|
||||
element.set("coeffs", ' '.join([str(self._coeffs[key])
|
||||
element.set("coeffs", ' '.join([str(self._coeffs.setdefault(key, 0.0))
|
||||
for key in self._coeff_keys]))
|
||||
|
||||
return element
|
||||
|
|
|
|||
|
|
@ -52,7 +52,7 @@ class Tally(object):
|
|||
List of nuclides to score results for
|
||||
scores : list of str
|
||||
List of defined scores, e.g. 'flux', 'fission', etc.
|
||||
estimator : {'analog', 'tracklength'}
|
||||
estimator : {'analog', 'tracklength', 'collision'}
|
||||
Type of estimator for the tally
|
||||
triggers : list of openmc.trigger.Trigger
|
||||
List of tally triggers
|
||||
|
|
@ -103,6 +103,9 @@ class Tally(object):
|
|||
self._with_batch_statistics = False
|
||||
self._derived = False
|
||||
|
||||
self._statepoint = None
|
||||
self._results_read = False
|
||||
|
||||
def __deepcopy__(self, memo):
|
||||
existing = memo.get(id(self))
|
||||
|
||||
|
|
@ -121,6 +124,8 @@ class Tally(object):
|
|||
clone._with_summary = self.with_summary
|
||||
clone._with_batch_statistics = self.with_batch_statistics
|
||||
clone._derived = self.derived
|
||||
clone._statepoint = self._statepoint
|
||||
clone._results_read = self._results_read
|
||||
|
||||
clone._filters = []
|
||||
for filter in self.filters:
|
||||
|
|
@ -265,24 +270,69 @@ class Tally(object):
|
|||
|
||||
@property
|
||||
def sum(self):
|
||||
if not self._statepoint:
|
||||
return None
|
||||
|
||||
if not self._results_read:
|
||||
# Extract Tally data from the file
|
||||
data = self._statepoint._f['tallies/tally {0}/results'.format(
|
||||
self.id)].value
|
||||
sum = data['sum']
|
||||
sum_sq = data['sum_sq']
|
||||
|
||||
# Define a routine to convert 0 to 1
|
||||
def nonzero(val):
|
||||
return 1 if not val else val
|
||||
|
||||
# Reshape the results arrays
|
||||
new_shape = (nonzero(self.num_filter_bins),
|
||||
nonzero(self.num_nuclides),
|
||||
nonzero(self.num_score_bins))
|
||||
|
||||
sum = np.reshape(sum, new_shape)
|
||||
sum_sq = np.reshape(sum_sq, new_shape)
|
||||
|
||||
# Set the data for this Tally
|
||||
self._sum = sum
|
||||
self._sum_sq = sum_sq
|
||||
|
||||
# Indicate that Tally results have been read
|
||||
self._results_read = True
|
||||
|
||||
return self._sum
|
||||
|
||||
@property
|
||||
def sum_sq(self):
|
||||
if not self._statepoint:
|
||||
return None
|
||||
|
||||
if not self._results_read:
|
||||
# Force reading of sum and sum_sq
|
||||
self.sum
|
||||
|
||||
return self._sum_sq
|
||||
|
||||
@property
|
||||
def mean(self):
|
||||
# Compute the mean if needed
|
||||
if self._mean is None:
|
||||
self.compute_mean()
|
||||
if not self._statepoint:
|
||||
return None
|
||||
|
||||
self._mean = self.sum / self.num_realizations
|
||||
return self._mean
|
||||
|
||||
@property
|
||||
def std_dev(self):
|
||||
# Compute the standard deviation if needed
|
||||
if self._std_dev is None:
|
||||
self.compute_std_dev()
|
||||
if not self._statepoint:
|
||||
return None
|
||||
|
||||
n = self.num_realizations
|
||||
nonzero = np.abs(self.mean) > 0
|
||||
self._std_dev = np.zeros_like(self.mean)
|
||||
self._std_dev[nonzero] = np.sqrt((self.sum_sq[nonzero]/n -
|
||||
self.mean[nonzero]**2)/(n - 1))
|
||||
self.with_batch_statistics = True
|
||||
return self._std_dev
|
||||
|
||||
@property
|
||||
|
|
@ -295,7 +345,8 @@ class Tally(object):
|
|||
|
||||
@estimator.setter
|
||||
def estimator(self, estimator):
|
||||
cv.check_value('estimator', estimator, ['analog', 'tracklength'])
|
||||
cv.check_value('estimator', estimator,
|
||||
['analog', 'tracklength', 'collision'])
|
||||
self._estimator = estimator
|
||||
|
||||
def add_trigger(self, trigger):
|
||||
|
|
@ -462,30 +513,6 @@ class Tally(object):
|
|||
|
||||
self._nuclides.remove(nuclide)
|
||||
|
||||
def compute_mean(self):
|
||||
"""Compute the sample mean for each bin in the tally"""
|
||||
|
||||
# Calculate sample mean
|
||||
self._mean = self.sum / self.num_realizations
|
||||
|
||||
def compute_std_dev(self, t_value=1.0):
|
||||
"""Compute the sample standard deviation for each bin in the tally
|
||||
|
||||
Parameters
|
||||
----------
|
||||
t_value : float, optional
|
||||
Student's t-value applied to the uncertainty. Defaults to 1.0,
|
||||
meaning the reported value is the sample standard deviation.
|
||||
|
||||
"""
|
||||
|
||||
# Calculate sample standard deviation
|
||||
self.compute_mean()
|
||||
self._std_dev = np.sqrt((self.sum_sq / self.num_realizations -
|
||||
self.mean**2) / (self.num_realizations - 1))
|
||||
self._std_dev *= t_value
|
||||
self.with_batch_statistics = True
|
||||
|
||||
def __repr__(self):
|
||||
string = 'Tally\n'
|
||||
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue