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

Merge branch 'develop' into arithmetic-units

Browse source
This commit is contained in:
Paul Romano 2017-10-24 13:27:39 -05:00
commit 5ca1d06b0c
45 changed files with 1308 additions and 868 deletions
Download patch file Download diff file Expand all files Collapse all files

10
docs/source/pythonapi/capi.rst
View file

@ -33,7 +33,9 @@ Classes
:nosignatures:
:template: myclass.rst
openmc.capi.CellView
openmc.capi.MaterialView
openmc.capi.NuclideView
openmc.capi.TallyView
openmc.capi.Cell
openmc.capi.EnergyFilter
openmc.capi.MaterialFilter
openmc.capi.Material
openmc.capi.Nuclide
openmc.capi.Tally

250
examples/jupyter/nuclear-data.ipynb
View file

File diff suppressed because one or more lines are too long

2
examples/python/pincell/build-xml.py
View file

@ -119,7 +119,7 @@ settings_file.export_to_xml()
###############################################################################
# Instantiate a tally mesh
mesh = openmc.Mesh(mesh_id=1)
mesh = openmc.Mesh()
mesh.type = 'regular'
mesh.dimension = [100, 100, 1]
mesh.lower_left = [-0.62992, -0.62992, -1.e50]

5
examples/xml/pincell/settings.xml
View file

@ -22,10 +22,11 @@
<!-- To assess convergence of the source distribution, we need to define the
bounds for a mesh over which the Shannon entropy should be
calculated. The extent in the z direction is made arbitrarily large. -->
<entropy>
<mesh id="1">
<lower_left>-0.39218 -0.39218 -1.e50</lower_left>
<upper_right>0.39218 0.39218 1.e50</upper_right>
<dimension>10 10 1</dimension>
</entropy>
</mesh>
<entropy_mesh>1</entropy_mesh>
</settings>

2
examples/xml/pincell/tallies.xml
View file

@ -1,7 +1,7 @@
<?xml version="1.0"?>
<tallies>
<mesh id="1" type="regular">
<mesh id="2" type="regular">
<dimension>100 100 1</dimension>
<lower_left>-0.62992 -0.62992 -1.e50</lower_left>
<upper_right>0.62992 0.62992 1.e50</upper_right>

12
openmc/arithmetic.py
View file

@ -108,9 +108,9 @@ class CrossNuclide(object):
Parameters
----------
left_nuclide : Nuclide or CrossNuclide
left_nuclide : openmc.Nuclide or CrossNuclide
The left nuclide in the outer product
right_nuclide : Nuclide or CrossNuclide
right_nuclide : openmc.Nuclide or CrossNuclide
The right nuclide in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
@ -118,9 +118,9 @@ class CrossNuclide(object):
Attributes
----------
left_nuclide : Nuclide or CrossNuclide
left_nuclide : openmc.Nuclide or CrossNuclide
The left nuclide in the outer product
right_nuclide : Nuclide or CrossNuclide
right_nuclide : openmc.Nuclide or CrossNuclide
The right nuclide in the outer product
binary_op : str
The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
@ -510,7 +510,7 @@ class AggregateNuclide(object):
Parameters
----------
nuclides : Iterable of str or Nuclide or CrossNuclide
nuclides : Iterable of str or openmc.Nuclide or CrossNuclide
The nuclides included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'avg', etc.) used
@ -518,7 +518,7 @@ class AggregateNuclide(object):
Attributes
----------
nuclides : Iterable of str or Nuclide or CrossNuclide
nuclides : Iterable of str or openmc.Nuclide or CrossNuclide
The nuclides included in the aggregation
aggregate_op : str
The tally aggregation operator (e.g., 'sum', 'avg', etc.) used

23
openmc/capi/__init__.py
View file

@ -13,6 +13,7 @@ objects in the :mod:`openmc.capi` subpackage, for example:
"""
from ctypes import CDLL
import os
import sys
from warnings import warn
@ -25,10 +26,21 @@ if sys.platform == 'darwin':
else:
_suffix = 'so'
# Open shared library
_filename = pkg_resources.resource_filename(
__name__, 'libopenmc.{}'.format(_suffix))
_dll = CDLL(_filename)
if os.environ.get('READTHEDOCS', None) != 'True':
# Open shared library
_filename = pkg_resources.resource_filename(
__name__, 'libopenmc.{}'.format(_suffix))
_dll = CDLL(_filename)
else:
# For documentation builds, we don't actually have the shared library
# available. Instead, we create a mock object so that when the modules
# within the openmc.capi package try to configure arguments and return
# values for symbols, no errors occur
try:
from unittest.mock import Mock
except ImportError:
from mock import Mock
_dll = Mock()
from .error import *
from .core import *
@ -38,3 +50,6 @@ from .cell import *
from .filter import *
from .tally import *
from .settings import settings
warn("The Python bindings to OpenMC's C API are still unstable "
"and may change substantially in future releases.", FutureWarning)

32
openmc/capi/cell.py
View file

@ -6,11 +6,11 @@ import numpy as np
from numpy.ctypeslib import as_array
from . import _dll
from .core import _View
from .core import _FortranObjectWithID
from .error import _error_handler
from .material import MaterialView
from .material import Material
__all__ = ['CellView', 'cells']
__all__ = ['Cell', 'cells']
# Cell functions
_dll.openmc_cell_get_id.argtypes = [c_int32, POINTER(c_int32)]
@ -31,11 +31,11 @@ _dll.openmc_get_cell_index.restype = c_int
_dll.openmc_get_cell_index.errcheck = _error_handler
class CellView(_View):
"""View of a cell.
class Cell(_FortranObjectWithID):
"""Cell stored internally.
This class exposes a cell that is stored internally in the OpenMC solver. To
obtain a view of a cell with a given ID, use the
This class exposes a cell that is stored internally in the OpenMC
library. To obtain a view of a cell with a given ID, use the
:data:`openmc.capi.nuclides` mapping.
Parameters
@ -53,7 +53,7 @@ class CellView(_View):
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
instance = super(Cell, self).__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
@ -72,9 +72,9 @@ class CellView(_View):
if fill_type.value == 1:
if n.value > 1:
return [MaterialView(i) for i in indices[:n.value]]
return [Material(index=i) for i in indices[:n.value]]
else:
return MaterialView(indices[0])
return Material(index=indices[0])
else:
raise NotImplementedError
@ -84,7 +84,7 @@ class CellView(_View):
n = len(fill)
indices = (c_int*n)(*(m._index for m in fill))
_dll.openmc_cell_set_fill(self._index, 1, 1, indices)
elif isinstance(fill, MaterialView):
elif isinstance(fill, Material):
materials = [fill]
indices = (c_int*1)(fill._index)
_dll.openmc_cell_set_fill(self._index, 1, 1, indices)
@ -108,12 +108,16 @@ class CellView(_View):
class _CellMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
_dll.openmc_get_cell_index(key, index)
return CellView(index.value)
try:
_dll.openmc_get_cell_index(key, index)
except (AllocationError, InvalidIDError) as e:
# __contains__ expects a KeyError to work correctly
raise KeyError(str(e))
return Cell(index.value)
def __iter__(self):
for i in range(len(self)):
yield CellView(i + 1).id
yield Cell(i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_cells').value

14
openmc/capi/core.py
View file

@ -172,9 +172,15 @@ class _DLLGlobal(object):
self.ctype.in_dll(_dll, self.name).value = value
class _View(object):
def __init__(self, index):
self._index = index
class _FortranObject(object):
def __repr__(self):
return "{}[{}]".format(type(self).__name__, self._index)
class _FortranObjectWithID(_FortranObject):
def __init__(self, uid=None, new=True, index=None):
# Creating the object has already been handled by __new__. In the
# initializer, all we do is make sure that the object returned has an ID
# assigned. If the array index of the object is out of bounds, an
# OutOfBoundsError will be raised here by virtue of referencing self.id
self.id

6
openmc/capi/error.py
View file

@ -39,9 +39,6 @@ class InvalidTypeError(Error):
"""Tried to perform an operation on the wrong type."""
_errmsg = (c_char*256).in_dll(_dll, 'openmc_err_msg')
def _error_handler(err, func, args):
"""Raise exception according to error code."""
@ -50,7 +47,8 @@ def _error_handler(err, func, args):
return c_int.in_dll(_dll, s).value
# Get error message set by OpenMC library
msg = _errmsg.value.decode()
errmsg = (c_char*256).in_dll(_dll, 'openmc_err_msg')
msg = errmsg.value.decode()
# Raise exception type corresponding to error code
if err == errcode('e_allocate'):

188
openmc/capi/filter.py
View file

@ -7,17 +7,17 @@ import numpy as np
from numpy.ctypeslib import as_array
from . import _dll
from .core import _View
from .error import _error_handler
from .material import MaterialView
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
from .material import Material
__all__ = ['FilterView', 'AzimuthalFilterView', 'CellFilterView',
'CellbornFilterView', 'CellfromFilterView', 'DistribcellFilterView',
'DelayedGroupFilterView', 'EnergyFilterView', 'EnergyoutFilterView',
'EnergyFunctionFilterView', 'MaterialFilterView', 'MeshFilterView',
'MuFilterView', 'PolarFilterView', 'SurfaceFilterView',
'UniverseFilterView', 'filters']
__all__ = ['Filter', 'AzimuthalFilter', 'CellFilter',
'CellbornFilter', 'CellfromFilter', 'DistribcellFilter',
'DelayedGroupFilter', 'EnergyFilter', 'EnergyoutFilter',
'EnergyFunctionFilter', 'MaterialFilter', 'MeshFilter',
'MuFilter', 'PolarFilter', 'SurfaceFilter',
'UniverseFilter', 'filters']
# Tally functions
_dll.openmc_energy_filter_get_bins.argtypes = [
@ -57,14 +57,43 @@ _dll.openmc_mesh_filter_set_mesh.restype = c_int
_dll.openmc_mesh_filter_set_mesh.errcheck = _error_handler
class FilterView(_View):
class Filter(_FortranObjectWithID):
__instances = WeakValueDictionary()
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
def __new__(cls, obj=None, uid=None, new=True, index=None):
mapping = filters
if index is None:
if new:
# Determine ID to assign
if uid is None:
try:
uid = max(mapping) + 1
except ValueError:
uid = 1
else:
if uid in mapping:
raise AllocationError('A filter with ID={} has already '
'been allocated.'.format(uid))
# Resize internal array
index = c_int32()
_dll.openmc_extend_filters(1, index, None)
# Set the filter type -- note that the filter_type attribute
# only exists on subclasses!
_dll.openmc_filter_set_type(index, cls.filter_type.encode())
index = index.value
else:
index = mapping[uid]._index
if index not in cls.__instances:
instance = super(Filter, cls).__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid
cls.__instances[index] = instance
return cls.__instances[index]
@property
def id(self):
@ -77,7 +106,14 @@ class FilterView(_View):
_dll.openmc_filter_set_id(self._index, filter_id)
class EnergyFilterView(FilterView):
class EnergyFilter(Filter):
filter_type = 'energy'
def __init__(self, bins=None, uid=None, new=True, index=None):
super(EnergyFilter, self).__init__(uid, new, index)
if bins is not None:
self.bins = bins
@property
def bins(self):
energies = POINTER(c_double)()
@ -95,45 +131,52 @@ class EnergyFilterView(FilterView):
self._index, len(energies), energies_p)
class EnergyoutFilterView(FilterView):
pass
class EnergyoutFilter(Filter):
filter_type = 'energyout'
class AzimuthalFilterView(FilterView):
pass
class AzimuthalFilter(Filter):
filter_type = 'azimuthal'
class CellFilterView(FilterView):
pass
class CellFilter(Filter):
filter_type = 'cell'
class CellbornFilterView(FilterView):
pass
class CellbornFilter(Filter):
filter_type = 'cellborn'
class CellfromFilterView(FilterView):
pass
class CellfromFilter(Filter):
filter_type = 'cellfrom'
class DelayedGroupFilterView(FilterView):
pass
class DelayedGroupFilter(Filter):
filter_type = 'delayedgroup'
class DistribcellFilterView(FilterView):
pass
class DistribcellFilter(Filter):
filter_type = 'distribcell'
class EnergyFunctionFilterView(FilterView):
pass
class EnergyFunctionFilter(Filter):
filter_type = 'energyfunction'
class MaterialFilterView(FilterView):
class MaterialFilter(Filter):
filter_type = 'material'
def __init__(self, bins=None, uid=None, new=True, index=None):
super(MaterialFilter, self).__init__(uid, new, index)
if bins is not None:
self.bins = bins
@property
def bins(self):
materials = POINTER(c_int32)()
n = c_int32()
_dll.openmc_material_filter_get_bins(self._index, materials, n)
return [MaterialView(materials[i]) for i in range(n.value)]
return [Material(index=materials[i]) for i in range(n.value)]
@bins.setter
def bins(self, materials):
@ -143,58 +186,43 @@ class MaterialFilterView(FilterView):
_dll.openmc_material_filter_set_bins(self._index, n, bins)
@classmethod
def new(cls, bins=None, filter_id=None):
# Determine filter ID to assign
if filter_id is None:
filter_id = max(filters) + 1
index = c_int32()
_dll.openmc_extend_filters(1, index, None)
_dll.openmc_filter_set_type(index, b'material')
f = cls(index.value)
f.id = filter_id
if bins is not None:
f.bins = bins
return f
class MeshFilter(Filter):
filter_type = 'mesh'
class MeshFilterView(FilterView):
pass
class MuFilter(Filter):
filter_type = 'mu'
class MuFilterView(FilterView):
pass
class PolarFilter(Filter):
filter_type = 'polar'
class PolarFilterView(FilterView):
pass
class SurfaceFilter(Filter):
filter_type = 'surface'
class SurfaceFilterView(FilterView):
pass
class UniverseFilter(Filter):
filter_type = 'universe'
class UniverseFilterView(FilterView):
pass
_filter_type_map = {
'azimuthal': AzimuthalFilterView,
'cell': CellFilterView,
'cellborn': CellbornFilterView,
'cellfrom': CellfromFilterView,
'delayedgroup': DelayedGroupFilterView,
'distribcell': DistribcellFilterView,
'energy': EnergyFilterView,
'energyout': EnergyoutFilterView,
'energyfunction': EnergyFunctionFilterView,
'material': MaterialFilterView,
'mesh': MeshFilterView,
'mu': MuFilterView,
'polar': PolarFilterView,
'surface': SurfaceFilterView,
'universe': UniverseFilterView,
_FILTER_TYPE_MAP = {
'azimuthal': AzimuthalFilter,
'cell': CellFilter,
'cellborn': CellbornFilter,
'cellfrom': CellfromFilter,
'delayedgroup': DelayedGroupFilter,
'distribcell': DistribcellFilter,
'energy': EnergyFilter,
'energyout': EnergyoutFilter,
'energyfunction': EnergyFunctionFilter,
'material': MaterialFilter,
'mesh': MeshFilter,
'mu': MuFilter,
'polar': PolarFilter,
'surface': SurfaceFilter,
'universe': UniverseFilter,
}
@ -202,18 +230,22 @@ def _get_filter(index):
filter_type = create_string_buffer(20)
_dll.openmc_filter_get_type(index, filter_type)
filter_type = filter_type.value.decode()
return _filter_type_map[filter_type](index)
return _FILTER_TYPE_MAP[filter_type](index=index)
class _FilterMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
_dll.openmc_get_filter_index(key, index)
try:
_dll.openmc_get_filter_index(key, index)
except (AllocationError, InvalidIDError) as e:
# __contains__ expects a KeyError to work correctly
raise KeyError(str(e))
return _get_filter(index.value)
def __iter__(self):
for i in range(len(self)):
yield TallyView(i + 1).id
yield _get_filter(i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_filters').value

75
openmc/capi/material.py
View file

@ -5,12 +5,12 @@ from weakref import WeakValueDictionary
import numpy as np
from numpy.ctypeslib import as_array
from . import _dll, NuclideView
from .core import _View
from .error import _error_handler
from . import _dll, Nuclide
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
__all__ = ['MaterialView', 'materials']
__all__ = ['Material', 'materials']
# Material functions
_dll.openmc_extend_materials.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)]
@ -43,16 +43,21 @@ _dll.openmc_material_set_id.restype = c_int
_dll.openmc_material_set_id.errcheck = _error_handler
class MaterialView(_View):
"""View of a material.
class Material(_FortranObjectWithID):
"""Material stored internally.
This class exposes a material that is stored internally in the OpenMC
solver. To obtain a view of a material with a given ID, use the
library. To obtain a view of a material with a given ID, use the
:data:`openmc.capi.materials` mapping.
Parameters
----------
index : int
uid : int or None
Unique ID of the tally
new : bool
When `index` is None, this argument controls whether a new object is
created or a view to an existing object is returned.
index : int or None
Index in the `materials` array.
Attributes
@ -67,6 +72,36 @@ class MaterialView(_View):
"""
__instances = WeakValueDictionary()
def __new__(cls, uid=None, new=True, index=None):
mapping = materials
if index is None:
if new:
# Determine ID to assign
if uid is None:
try:
uid = max(mapping) + 1
except ValueError:
uid = 1
else:
if uid in mapping:
raise AllocationError('A material with ID={} has already '
'been allocated.'.format(uid))
index = c_int32()
_dll.openmc_extend_materials(1, index, None)
index = index.value
else:
index = mapping[uid]._index
if index not in cls.__instances:
instance = super(Material, cls).__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid
cls.__instances[index] = instance
return cls.__instances[index]
@property
def id(self):
mat_id = c_int32()
@ -106,7 +141,7 @@ class MaterialView(_View):
_dll.openmc_material_get_densities(self._index, nuclides, densities, n)
# Convert to appropriate types and return
nuclide_list = [NuclideView(nuclides[i]).name for i in range(n.value)]
nuclide_list = [Nuclide(nuclides[i]).name for i in range(n.value)]
density_array = as_array(densities, (n.value,))
return nuclide_list, density_array
@ -123,18 +158,6 @@ class MaterialView(_View):
"""
_dll.openmc_material_add_nuclide(self._index, name.encode(), density)
@classmethod
def new(cls, material_id=None):
# Determine ID to assign
if material_id is None:
material_id = max(materials) + 1
index = c_int32()
_dll.openmc_extend_materials(1, index, None)
mat = cls(index.value)
mat.id = material_id
return mat
def set_density(self, density):
"""Set density of a material.
@ -171,12 +194,16 @@ class MaterialView(_View):
class _MaterialMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
_dll.openmc_get_material_index(key, index)
return MaterialView(index.value)
try:
_dll.openmc_get_material_index(key, index)
except (AllocationError, InvalidIDError) as e:
# __contains__ expects a KeyError to work correctly
raise KeyError(str(e))
return Material(index=index.value)
def __iter__(self):
for i in range(len(self)):
yield MaterialView(i + 1).id
yield Material(index=i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_materials').value

25
openmc/capi/nuclide.py
View file

@ -6,11 +6,11 @@ import numpy as np
from numpy.ctypeslib import as_array
from . import _dll
from .core import _View
from .error import _error_handler
from .core import _FortranObject
from .error import _error_handler, DataError, AllocationError
__all__ = ['NuclideView', 'nuclides', 'load_nuclide']
__all__ = ['Nuclide', 'nuclides', 'load_nuclide']
# Nuclide functions
_dll.openmc_get_nuclide_index.argtypes = [c_char_p, POINTER(c_int)]
@ -36,8 +36,8 @@ def load_nuclide(name):
_dll.openmc_load_nuclide(name.encode())
class NuclideView(_View):
"""View of a nuclide.
class Nuclide(_FortranObject):
"""Nuclide stored internally.
This class exposes a nuclide that is stored internally in the OpenMC
solver. To obtain a view of a nuclide with a given name, use the
@ -58,10 +58,13 @@ class NuclideView(_View):
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
instance = super(Nuclide, cls).__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
def __init__(self, index):
self._index = index
@property
def name(self):
name = c_char_p()
@ -78,12 +81,16 @@ class _NuclideMapping(Mapping):
"""Provide mapping from nuclide name to index in nuclides array."""
def __getitem__(self, key):
index = c_int()
_dll.openmc_get_nuclide_index(key.encode(), index)
return NuclideView(index.value)
try:
_dll.openmc_get_nuclide_index(key.encode(), index)
except (DataError, AllocationError) as e:
# __contains__ expects a KeyError to work correctly
raise KeyError(str(e))
return Nuclide(index.value)
def __iter__(self):
for i in range(len(self)):
yield NuclideView(i + 1).name
yield Nuclide(i + 1).name
def __len__(self):
return c_int.in_dll(_dll, 'n_nuclides').value

12
openmc/capi/settings.py
View file

@ -10,6 +10,10 @@ _RUN_MODES = {1: 'fixed source',
4: 'particle restart',
5: 'volume'}
_dll.openmc_set_seed.argtypes = [c_int64]
_dll.openmc_set_seed.restype = c_int
_dll.openmc_set_seed.errcheck = _error_handler
class _Settings(object):
# Attributes that are accessed through a descriptor
@ -37,5 +41,13 @@ class _Settings(object):
else:
raise ValueError('Invalid run mode: {}'.format(mode))
@property
def seed(self):
return c_int64.in_dll(_dll, 'seed').value
@seed.setter
def seed(self, seed):
_dll.openmc_set_seed(seed)
settings = _Settings()

77
openmc/capi/tally.py
View file

@ -4,13 +4,13 @@ from weakref import WeakValueDictionary
from numpy.ctypeslib import as_array
from . import _dll, NuclideView
from .core import _View
from .error import _error_handler
from . import _dll, Nuclide
from .core import _FortranObjectWithID
from .error import _error_handler, AllocationError, InvalidIDError
from .filter import _get_filter
__all__ = ['TallyView', 'tallies']
__all__ = ['Tally', 'tallies']
# Tally functions
_dll.openmc_get_tally_index.argtypes = [c_int32, POINTER(c_int32)]
@ -51,24 +51,29 @@ _dll.openmc_tally_set_type.restype = c_int
_dll.openmc_tally_set_type.errcheck = _error_handler
class TallyView(_View):
"""View of a tally.
class Tally(_FortranObjectWithID):
"""Tally stored internally.
This class exposes a tally that is stored internally in the OpenMC
solver. To obtain a view of a tally with a given ID, use the
library. To obtain a view of a tally with a given ID, use the
:data:`openmc.capi.tallies` mapping.
Parameters
----------
index : int
Index in the `tallies` array.
uid : int or None
Unique ID of the tally
new : bool
When `index` is None, this argument controls whether a new object is
created or a view of an existing object is returned.
index : int or None
Index in the `tallies` array.
Attributes
----------
id : int
ID of the tally
filters : list
List of views to tally filters
List of tally filters
nuclides : list of str
List of nuclides to score results for
results : numpy.ndarray
@ -77,11 +82,36 @@ class TallyView(_View):
"""
__instances = WeakValueDictionary()
def __new__(cls, *args):
if args not in cls.__instances:
instance = super().__new__(cls)
cls.__instances[args] = instance
return cls.__instances[args]
def __new__(cls, uid=None, new=True, index=None):
mapping = tallies
if index is None:
if new:
# Determine ID to assign
if uid is None:
try:
uid = max(mapping) + 1
except ValueError:
uid = 1
else:
if uid in mapping:
raise AllocationError('A tally with ID={} has already '
'been allocated.'.format(uid))
index = c_int32()
_dll.openmc_extend_tallies(1, index, None)
_dll.openmc_tally_set_type(index, b'generic')
index = index.value
else:
index = mapping[uid]._index
if index not in cls.__instances:
instance = super(Tally, cls).__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid
cls.__instances[index] = instance
return cls.__instances[index]
@property
def id(self):
@ -105,7 +135,7 @@ class TallyView(_View):
nucs = POINTER(c_int)()
n = c_int()
_dll.openmc_tally_get_nuclides(self._index, nucs, n)
return [NuclideView(nucs[i]).name if nucs[i] > 0 else 'total'
return [Nuclide(nucs[i]).name if nucs[i] > 0 else 'total'
for i in range(n.value)]
@property
@ -143,7 +173,10 @@ class TallyView(_View):
def new(cls, tally_id=None):
# Determine ID to assign
if tally_id is None:
tally_id = max(tallies) + 1
try:
tally_id = max(tallies) + 1
except ValueError:
tally_id = 1
index = c_int32()
_dll.openmc_extend_tallies(1, index, None)
@ -156,12 +189,16 @@ class TallyView(_View):
class _TallyMapping(Mapping):
def __getitem__(self, key):
index = c_int32()
_dll.openmc_get_tally_index(key, index)
return TallyView(index.value)
try:
_dll.openmc_get_tally_index(key, index)
except (AllocationError, InvalidIDError) as e:
# __contains__ expects a KeyError to work correctly
raise KeyError(str(e))
return Tally(index=index.value)
def __iter__(self):
for i in range(len(self)):
yield TallyView(i + 1).id
yield Tally(index=i + 1).id
def __len__(self):
return c_int32.in_dll(_dll, 'n_tallies').value

14
openmc/element.py
View file

@ -124,6 +124,15 @@ class Element(object):
is a tuple consisting of an openmc.Nuclide instance and the natural
abundance of the isotope.
Notes
-----
When the `enrichment` argument is specified, a correlation from
`ORNL/CSD/TM-244 <https://doi.org/10.2172/5561567>`_ is used to
calculate the weight fractions of U234, U235, U236, and U238. Namely,
the weight fraction of U234 and U236 are taken to be 0.89% and 0.46%,
respectively, of the U235 weight fraction. The remainder of the isotopic
weight is assigned to U238.
"""
# Get the nuclides present in nature
@ -217,9 +226,10 @@ class Element(object):
if enrichment is not None:
# Calculate the mass fractions of isotopes
abundances['U234'] = 0.008 * enrichment
abundances['U234'] = 0.0089 * enrichment
abundances['U235'] = enrichment
abundances['U238'] = 100.0 - 1.008 * enrichment
abundances['U236'] = 0.0046 * enrichment
abundances['U238'] = 100.0 - 1.0135 * enrichment
# Convert the mass fractions to mole fractions
for nuclide in abundances.keys():

6
openmc/mgxs/library.py
View file

@ -606,7 +606,7 @@ class Library(object):
Parameters
----------
domain : Material or Cell or Universe or Integral
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh or Integral
The material, cell, or universe object of interest (or its ID)
mgxs_type : {'total', 'transport', 'nu-transport', 'absorption', 'capture', 'fission', 'nu-fission', 'kappa-fission', 'scatter', 'nu-scatter', 'scatter matrix', 'nu-scatter matrix', 'multiplicity matrix', 'nu-fission matrix', chi', 'chi-prompt', 'inverse-velocity', 'prompt-nu-fission', 'prompt-nu-fission matrix', 'delayed-nu-fission', 'delayed-nu-fission matrix', 'chi-delayed', 'beta'}
The type of multi-group cross section object to return
@ -668,7 +668,7 @@ class Library(object):
Returns
-------
Library
openmc.mgxs.Library
A new multi-group cross section library condensed to the group
structure of interest
@ -880,7 +880,7 @@ class Library(object):
Returns
-------
Library
openmc.mgxs.Library
A Library object loaded from the pickle binary file
See also

14
openmc/mgxs/mdgxs.py
View file

@ -71,7 +71,7 @@ class MDGXS(MGXS):
Reaction type (e.g., 'chi-delayed', 'beta', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -948,7 +948,7 @@ class ChiDelayed(MDGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -1462,7 +1462,7 @@ class DelayedNuFissionXS(MDGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -1598,7 +1598,7 @@ class Beta(MDGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -1782,7 +1782,7 @@ class DecayRate(MDGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -1959,7 +1959,7 @@ class MatrixMDGXS(MDGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -2551,7 +2551,7 @@ class DelayedNuFissionMatrixXS(MatrixMDGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization

32
openmc/mgxs/mgxs.py
View file

@ -154,7 +154,7 @@ class MGXS(object):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -2018,7 +2018,7 @@ class MatrixMGXS(MGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -2512,7 +2512,7 @@ class TotalXS(MGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -2649,7 +2649,7 @@ class TransportXS(MGXS):
If True, the cross section data will include neutron multiplication
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -2858,7 +2858,7 @@ class AbsorptionXS(MGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -2986,7 +2986,7 @@ class CaptureXS(MGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -3140,7 +3140,7 @@ class FissionXS(MGXS):
If true, computes cross sections which only includes prompt neutrons
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -3309,7 +3309,7 @@ class KappaFissionXS(MGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -3441,7 +3441,7 @@ class ScatterXS(MGXS):
If True, the cross section data will include neutron multiplication
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -3658,7 +3658,7 @@ class ScatterMatrixXS(MatrixMGXS):
If True, the cross section data will include neutron multiplication
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -4020,7 +4020,7 @@ class ScatterMatrixXS(MatrixMGXS):
# Override the nuclides for tally arithmetic
correction.nuclides = scatter_p1.nuclides
self._xs_tally -= correction
self._compute_xs()
return self._xs_tally
@ -4749,7 +4749,7 @@ class MultiplicityMatrixXS(MatrixMGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -4915,7 +4915,7 @@ class ScatterProbabilityMatrix(MatrixMGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -5088,7 +5088,7 @@ class NuFissionMatrixXS(MatrixMGXS):
If true, computes cross sections which only includes prompt neutrons
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -5245,7 +5245,7 @@ class Chi(MGXS):
If true, computes cross sections which only includes prompt neutrons
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization
@ -5823,7 +5823,7 @@ class InverseVelocity(MGXS):
Reaction type (e.g., 'total', 'nu-fission', etc.)
by_nuclide : bool
If true, computes cross sections for each nuclide in domain
domain : Material or Cell or Universe or Mesh
domain : openmc.Material or openmc.Cell or openmc.Universe or openmc.Mesh
Domain for spatial homogenization
domain_type : {'material', 'cell', 'distribcell', 'universe', 'mesh'}
Domain type for spatial homogenization

2
openmc/tallies.py
View file

@ -502,7 +502,7 @@ class Tally(IDManagerMixin):
Parameters
----------
nuclide : str, Nuclide, CrossNuclide or AggregateNuclide
nuclide : str, openmc.Nuclide, CrossNuclide or AggregateNuclide
Nuclide to add to the tally. The nuclide should be a Nuclide object
when a user is adding nuclides to a Tally for input file generation.
The nuclide is a str when a Tally is created from a StatePoint file

8
src/api.F90
View file

@ -17,7 +17,7 @@ module openmc_api
use initialize, only: openmc_init
use particle_header, only: Particle
use plot, only: openmc_plot_geometry
use random_lcg, only: seed, initialize_prng
use random_lcg, only: seed, openmc_set_seed
use settings
use simulation_header
use tally_header
@ -51,6 +51,7 @@ module openmc_api
public :: openmc_get_cell_index
public :: openmc_get_keff
public :: openmc_get_filter_index
public :: openmc_get_filter_next_id
public :: openmc_get_material_index
public :: openmc_get_nuclide_index
public :: openmc_get_tally_index
@ -212,6 +213,8 @@ contains
!===============================================================================
subroutine openmc_hard_reset() bind(C)
integer :: err
! Reset all tallies and timers
call openmc_reset()
@ -220,8 +223,7 @@ contains
total_gen = 0
! Reset the random number generator state
seed = 1_8
call initialize_prng()
err = openmc_set_seed(1_8)
end subroutine openmc_hard_reset
!===============================================================================

19
src/cmfd_input.F90
View file

@ -261,7 +261,8 @@ contains
integer :: i_filt_start, i_filt_end
integer(C_INT32_T), allocatable :: filter_indices(:)
integer(C_INT) :: err
integer :: i_filt ! index in filters array
integer :: i_filt ! index in filters array
integer :: filt_id
integer :: iarray3(3) ! temp integer array
real(8) :: rarray3(3) ! temp double array
real(C_DOUBLE), allocatable :: energies(:)
@ -369,7 +370,7 @@ contains
m % volume_frac = ONE/real(product(m % dimension),8)
! Add mesh to dictionary
call mesh_dict % add_key(m % id, i_start)
call mesh_dict % set(m % id, i_start)
! Determine number of filters
energy_filters = check_for_node(node_mesh, "energy")
@ -381,14 +382,16 @@ contains
! Set up mesh filter
i_filt = i_filt_start
err = openmc_filter_set_type(i_filt, C_CHAR_'mesh' // C_NULL_CHAR)
err = openmc_filter_set_id(i_filt, i_filt)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_mesh_filter_set_mesh(i_filt, i_start)
if (energy_filters) then
! Read and set incoming energy mesh filter
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'energy' // C_NULL_CHAR)
err = openmc_filter_set_id(i_filt, i_filt)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
! Get energies and set bins
ng = node_word_count(node_mesh, "energy")
@ -399,7 +402,8 @@ contains
! Read and set outgoing energy mesh filter
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'energyout' // C_NULL_CHAR)
err = openmc_filter_set_id(i_filt, i_filt)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_energy_filter_set_bins(i_filt, ng, energies)
end if
@ -407,7 +411,8 @@ contains
! tallies use this filter and we need to change the dimension
i_filt = i_filt + 1
err = openmc_filter_set_type(i_filt, C_CHAR_'mesh' // C_NULL_CHAR)
err = openmc_filter_set_id(i_filt, i_filt)
call openmc_get_filter_next_id(filt_id)
err = openmc_filter_set_id(i_filt, filt_id)
err = openmc_mesh_filter_set_mesh(i_filt, i_start)
! We need to increase the dimension by one since we also need
@ -432,7 +437,7 @@ contains
end if
filt % current = .true.
! Add filter to dictionary
call filter_dict % add_key(filt % id, i_filt)
call filter_dict % set(filt % id, i_filt)
end select
! Initialize filters

767
src/dict_header.F90
View file

@ -7,438 +7,535 @@ module dict_header
! pairs. This data structure is used to provide lookup features, e.g. cells and
! surfaces by name.
!
! The original version was roughly based on capabilities in the 'flibs' open
! source package. However, it was rewritten from scratch so that it could be
! used stand-alone without relying on the implementation of lists. As with
! lists, it was considered writing a single dictionary used unlimited
! polymorphism, but again compiler support is spotty and doesn't always prevent
! duplication of code.
! The implementation is based on Algorithm D from Knuth Vol. 3 Sec. 6.4 (open
! addressing with double hashing). Hash table sizes M are chosen such that M
! and M - 2 are twin primes, which helps reduce clustering. The sequence of
! twin primes used for the table sizes comes from
! https://github.com/anholt/hash_table/blob/master/hash_table.c. These values
! were selected so that the table would grow by approximately a factor of two
! each time the maximum load factor is exceeded. An upper limit is placed on
! the load factor to prevent exponential performance degradation as the number
! of entries approaches the number of buckets.
!===============================================================================
implicit none
integer, parameter, private :: HASH_SIZE = 4993
integer, parameter, private :: HASH_MULTIPLIER = 31
integer, parameter, private :: DICT_NULL = -huge(0)
integer, parameter :: DICT_KEY_LENGTH = 255
integer, parameter :: EMPTY = -huge(0)
integer, parameter, private :: DELETED = -huge(0) + 1
integer, parameter, private :: KEY_CHAR_LENGTH = 255
integer, parameter, private :: TABLE_SIZES(30) = &
[5, 7, 13, 19, 43, 73, 151, 283, 571, 1153, 2269, 4519, 9013, 18043, &
36109, 72091, 144409, 288361, 576883, 1153459, 2307163, 4613893, &
9227641, 18455029, 36911011, 73819861, 147639589, 295279081, &
590559793, 1181116273]
real(8), parameter, private :: MAX_LOAD_FACTOR = 0.65
!===============================================================================
! ELEMKEYVALUE* contains (key,value) pairs and a pointer to the next (key,value)
! pair
! DICTENTRY* contains (key,value) pairs.
!===============================================================================
type ElemKeyValueCI
type(ElemKeyValueCI), pointer :: next => null()
character(len=DICT_KEY_LENGTH) :: key
integer :: value
end type ElemKeyValueCI
type DictEntryII
integer :: key = EMPTY
integer :: value = EMPTY
end type DictEntryII
type ElemKeyValueII
type(ElemKeyValueII), pointer :: next => null()
integer :: key
integer :: value
end type ElemKeyValueII
type DictEntryCI
character(len=KEY_CHAR_LENGTH) :: key
integer :: value = EMPTY
end type DictEntryCI
!===============================================================================
! HASHLIST* types contain a single pointer to a linked list of (key,value)
! pairs. This type is necesssary so that the Dict types can be dynamically
! allocated.
! BUCKET* contains an allocatable DictEntry object for storing a (key,value)
! pair and the hash value for fast comparisons. Integer (key,value) pairs are
! stored directly in the hash table since their memory requirement is small.
!===============================================================================
type, private :: HashListCI
type(ElemKeyValueCI), pointer :: list => null()
end type HashListCI
type, private :: HashListII
type(ElemKeyValueII), pointer :: list => null()
end type HashListII
type, private :: BucketCI
type(DictEntryCI), allocatable :: entry
integer :: hash = EMPTY
end type BucketCI
!===============================================================================
! DICT* is a dictionary of (key,value) pairs with convenience methods as
! type-bound procedures. DictCharInt has character(*) keys and integer values,
! and DictIntInt has integer keys and values.
! type-bound procedures. DictIntInt has integer keys and values, and
! DictCharInt has character(*) keys and integer values.
!===============================================================================
type, public :: DictCharInt
private
type(HashListCI), pointer :: table(:) => null()
contains
procedure :: add_key => dict_add_key_ci
procedure :: get_key => dict_get_key_ci
procedure :: has_key => dict_has_key_ci
procedure :: keys => dict_keys_ci
procedure :: clear => dict_clear_ci
procedure, private :: get_elem => dict_get_elem_ci
end type DictCharInt
type, public :: DictIntInt
private
type(HashListII), pointer :: table(:) => null()
integer, private :: entries = 0
integer, private :: capacity = 0
type(DictEntryII), allocatable, private :: table(:)
contains
procedure :: add_key => dict_add_key_ii
procedure :: get_key => dict_get_key_ii
procedure :: has_key => dict_has_key_ii
procedure :: keys => dict_keys_ii
procedure :: clear => dict_clear_ii
procedure, private :: get_elem => dict_get_elem_ii
procedure :: set => set_ii
procedure :: get => get_ii
procedure :: has => has_ii
procedure :: remove => remove_ii
procedure :: next_entry => next_entry_ii
procedure :: clear => clear_ii
procedure :: size => size_ii
procedure, private :: get_entry => get_entry_ii
procedure, private :: resize => resize_ii
end type DictIntInt
type, public :: DictCharInt
integer, private :: entries = 0
integer, private :: capacity = 0
type(BucketCI), allocatable, private :: table(:)
contains
procedure :: set => set_ci
procedure :: get => get_ci
procedure :: has => has_ci
procedure :: remove => remove_ci
procedure :: next_entry => next_entry_ci
procedure :: clear => clear_ci
procedure :: size => size_ci
procedure, private :: get_entry => get_entry_ci
procedure, private :: resize => resize_ci
end type DictCharInt
contains
!===============================================================================
! DICT_ADD_KEY adds a (key,value) entry to a dictionary. If the key is already
! in the dictionary, the value is replaced by the new specified value.
! GET_ENTRY returns the index of the (key,value) pair in the table for a given
! key. This method is private.
!===============================================================================
subroutine dict_add_key_ci(this, key, value)
function get_entry_ii(this, key) result(i)
class(DictCharInt) :: this
character(*), intent(in) :: key
integer, intent(in) :: value
class(DictIntInt) :: this
integer, intent(in) :: key
integer :: i
integer :: hash
type(ElemKeyValueCI), pointer :: elem => null()
type(ElemKeyValueCI), pointer :: new_elem => null()
integer :: c
elem => this % get_elem(key)
if (associated(elem)) then
elem % value = value
else
! Get hash
hash = dict_hash_key_ci(key)
! Create new element
allocate(new_elem)
new_elem % key = key
new_elem % value = value
! Add element to front of list
new_elem % next => this % table(hash) % list
this % table(hash) % list => new_elem
if (.not. allocated(this % table)) then
allocate(this % table(TABLE_SIZES(1)))
this % capacity = TABLE_SIZES(1)
end if
end subroutine dict_add_key_ci
hash = hash_ii(key)
i = 1 + mod(hash, this % capacity)
c = 2 + mod(hash, this % capacity - 2)
subroutine dict_add_key_ii(this, key, value)
do
if (this % table(i) % key == key .or. &
this % table(i) % key == EMPTY) exit
class(DictIntInt) :: this
i = 1 + mod(i + c - 1, this % capacity)
end do
end function get_entry_ii
function get_entry_ci(this, key) result(i)
class(DictCharInt) :: this
character(*), intent(in) :: key
integer :: i
integer :: hash
integer :: c
if (.not. allocated(this % table)) then
allocate(this % table(TABLE_SIZES(1)))
this % capacity = TABLE_SIZES(1)
end if
hash = hash_ci(key)
i = 1 + mod(hash, this % capacity)
c = 2 + mod(hash, this % capacity - 2)
do
if (this % table(i) % hash == hash) then
if (allocated(this % table(i) % entry)) then
if (this % table(i) % entry % key == key) exit
end if
end if
if (this % table(i) % hash == EMPTY) exit
i = 1 + mod(i + c - 1, this % capacity)
end do
end function get_entry_ci
!===============================================================================
! RESIZE allocates a new hash table to accomodate the number of entries and
! reinserts all of the entries into the new table. This method is private.
!===============================================================================
subroutine resize_ii(this)
class(DictIntInt) :: this
type(DictEntryII), allocatable :: table(:)
integer :: new_size
integer :: i
do i = 1, size(TABLE_SIZES)
if (TABLE_SIZES(i) > this % capacity) exit
end do
new_size = TABLE_SIZES(i)
call move_alloc(this % table, table)
allocate(this % table(new_size))
this % capacity = new_size
this % entries = 0
! Rehash each entry into the new table
do i = 1, size(table)
if (table(i) % key /= EMPTY .and. table(i) % key /= DELETED) then
call this % set(table(i) % key, table(i) % value)
end if
end do
deallocate(table)
end subroutine resize_ii
subroutine resize_ci(this)
class(DictCharInt) :: this
type(BucketCI), allocatable :: table(:)
integer :: new_size
integer :: i
do i = 1, size(TABLE_SIZES)
if (TABLE_SIZES(i) > this % capacity) exit
end do
new_size = TABLE_SIZES(i)
call move_alloc(this % table, table)
allocate(this % table(new_size))
this % capacity = new_size
this % entries = 0
! Rehash each entry into the new table
do i = 1, size(table)
if (table(i) % hash /= EMPTY .and. table(i) % hash /= DELETED) then
call this % set(table(i) % entry % key, table(i) % entry % value)
end if
end do
deallocate(table)
end subroutine resize_ci
!===============================================================================
! SET adds a (key,value) entry to a dictionary. If the key is already in the
! dictionary, the value is replaced by the new specified value.
!===============================================================================
subroutine set_ii(this, key, value)
class(DictIntInt) :: this
integer, intent(in) :: key
integer, intent(in) :: value
integer :: hash
type(ElemKeyValueII), pointer :: elem => null()
type(ElemKeyValueII), pointer :: new_elem => null()
integer :: i
integer :: c
elem => this % get_elem(key)
if (associated(elem)) then
elem % value = value
else
! Get hash
hash = dict_hash_key_ii(key)
! Create new element
allocate(new_elem)
new_elem % key = key
new_elem % value = value
! Add element to front of list
new_elem % next => this % table(hash) % list
this % table(hash) % list => new_elem
if (.not. allocated(this % table)) then
allocate(this % table(TABLE_SIZES(1)))
this % capacity = TABLE_SIZES(1)
else if (real(this % entries + 1, 8) / this % capacity > MAX_LOAD_FACTOR) then
call this % resize()
end if
end subroutine dict_add_key_ii
hash = hash_ii(key)
i = 1 + mod(hash, this % capacity)
c = 2 + mod(hash, this % capacity - 2)
!===============================================================================
! DICT_GET_ELEM returns a pointer to the (key,value) pair for a given key. This
! method is private.
!===============================================================================
do
if (this % table(i) % key == EMPTY .or. &
this % table(i) % key == DELETED) then
this % table(i) % key = key
this % table(i) % value = value
this % entries = this % entries + 1
exit
else if (this % table(i) % key == key) then
this % table(i) % value = value
exit
end if
function dict_get_elem_ci(this, key) result(elem)
class(DictCharInt) :: this
character(*), intent(in) :: key
type(ElemKeyValueCI), pointer :: elem
integer :: hash
! Check for dictionary not being allocated
if (.not. associated(this % table)) then
allocate(this % table(HASH_SIZE))
end if
hash = dict_hash_key_ci(key)
elem => this % table(hash) % list
do while (associated(elem))
if (elem % key == key) exit
elem => elem % next
i = 1 + mod(i + c - 1, this % capacity)
end do
end function dict_get_elem_ci
end subroutine set_ii
function dict_get_elem_ii(this, key) result(elem)
subroutine set_ci(this, key, value)
class(DictIntInt) :: this
integer, intent(in) :: key
type(ElemKeyValueII), pointer :: elem
integer :: hash
! Check for dictionary not being allocated
if (.not. associated(this % table)) then
allocate(this % table(HASH_SIZE))
end if
hash = dict_hash_key_ii(key)
elem => this % table(hash) % list
do while (associated(elem))
if (elem % key == key) exit
elem => elem % next
end do
end function dict_get_elem_ii
!===============================================================================
! DICT_GET_KEY returns the value matching a given key. If the dictionary does
! not contain the key, the value DICT_NULL is returned.
!===============================================================================
function dict_get_key_ci(this, key) result(value)
class(DictCharInt) :: this
class(DictCharInt) :: this
character(*), intent(in) :: key
integer :: value
integer, intent(in) :: value
type(ElemKeyValueCI), pointer :: elem
integer :: hash
integer :: i
integer :: c
elem => this % get_elem(key)
if (associated(elem)) then
value = elem % value
else
value = DICT_NULL
if (.not. allocated(this % table)) then
allocate(this % table(TABLE_SIZES(1)))
this % capacity = TABLE_SIZES(1)
else if (real(this % entries + 1, 8) / this % capacity > MAX_LOAD_FACTOR) then
call this % resize()
end if
end function dict_get_key_ci
hash = hash_ci(key)
i = 1 + mod(hash, this % capacity)
c = 2 + mod(hash, this % capacity - 2)
function dict_get_key_ii(this, key) result(value)
do
if (this % table(i) % hash == EMPTY .or. &
this % table(i) % hash == DELETED) then
if (.not. allocated(this % table(i) % entry)) then
allocate(this % table(i) % entry)
end if
this % table(i) % hash = hash
this % table(i) % entry % key = key
this % table(i) % entry % value = value
this % entries = this % entries + 1
exit
else if (this % table(i) % hash == hash .and. &
this % table(i) % entry % key == key) then
this % table(i) % entry % value = value
exit
end if
i = 1 + mod(i + c - 1, this % capacity)
end do
end subroutine set_ci
!===============================================================================
! GET returns the value matching a given key. If the dictionary does not contain
! the key, the value EMPTY is returned.
!===============================================================================
function get_ii(this, key) result(value)
class(DictIntInt) :: this
integer, intent(in) :: key
integer :: value
type(ElemKeyValueII), pointer :: elem
integer :: i
elem => this % get_elem(key)
i = this % get_entry(key)
value = this % table(i) % value
if (associated(elem)) then
value = elem % value
end function get_ii
function get_ci(this, key) result(value)
class(DictCharInt) :: this
character(*), intent(in) :: key
integer :: value
integer :: i
i = this % get_entry(key)
if (allocated(this % table(i) % entry)) then
value = this % table(i) % entry % value
else
value = DICT_NULL
value = EMPTY
end if
end function dict_get_key_ii
end function get_ci
!===============================================================================
! DICT_HAS_KEY determines whether a dictionary has a (key,value) pair with a
! given key.
! HAS determines whether a dictionary has a (key,value) pair with a given key.
!===============================================================================
function dict_has_key_ci(this, key) result(has)
class(DictCharInt) :: this
character(*), intent(in) :: key
logical :: has
type(ElemKeyValueCI), pointer :: elem
elem => this % get_elem(key)
has = associated(elem)
end function dict_has_key_ci
function dict_has_key_ii(this, key) result(has)
function has_ii(this, key) result(has)
class(DictIntInt) :: this
integer, intent(in) :: key
logical :: has
type(ElemKeyValueII), pointer :: elem
integer :: i
elem => this % get_elem(key)
has = associated(elem)
i = this % get_entry(key)
has = (this % table(i) % key /= EMPTY)
end function dict_has_key_ii
end function has_ii
!===============================================================================
! DICT_HASH_KEY returns the hash value for a given key
!===============================================================================
function dict_hash_key_ci(key) result(val)
function has_ci(this, key) result(has)
class(DictCharInt) :: this
character(*), intent(in) :: key
integer :: val
logical :: has
integer :: i
val = 0
i = this % get_entry(key)
has = (this % table(i) % hash /= EMPTY)
do i = 1, len_trim(key)
val = HASH_MULTIPLIER * val + ichar(key(i:i))
end do
end function has_ci
! Added the absolute val on val-1 since the sum in the do loop is
! susceptible to integer overflow
val = 1 + mod(abs(val-1), HASH_SIZE)
!===============================================================================
! REMOVE deletes a (key,value) entry from a dictionary.
!===============================================================================
end function dict_hash_key_ci
function dict_hash_key_ii(key) result(val)
subroutine remove_ii(this, key)
class(DictIntInt) :: this
integer, intent(in) :: key
integer :: val
val = 0
! Added the absolute val on val-1 since the sum in the do loop is
! susceptible to integer overflow
val = 1 + mod(abs(key-1), HASH_SIZE)
end function dict_hash_key_ii
!===============================================================================
! DICT_KEYS returns a pointer to a linked list of all (key,value) pairs
!===============================================================================
function dict_keys_ci(this) result(keys)
class(DictCharInt) :: this
type(ElemKeyValueCI), pointer :: keys
integer :: i
type(ElemKeyValueCI), pointer :: current => null()
type(ElemKeyValueCI), pointer :: elem => null()
keys => null()
do i = 1, size(this % table)
! Get pointer to start of bucket i
elem => this % table(i) % list
do while (associated(elem))
! Allocate (key,value) pair
if (.not. associated(keys)) then
allocate(keys)
current => keys
else
allocate(current % next)
current => current % next
end if
! Copy (key,value) pair
current % key = elem % key
current % value = elem % value
! Move to next element in bucket i
elem => elem % next
end do
end do
end function dict_keys_ci
function dict_keys_ii(this) result(keys)
class(DictIntInt) :: this
type(ElemKeyValueII), pointer :: keys
integer :: i
type(ElemKeyValueII), pointer :: current => null()
type(ElemKeyValueII), pointer :: elem => null()
keys => null()
do i = 1, size(this % table)
! Get pointer to start of bucket i
elem => this % table(i) % list
do while (associated(elem))
! Allocate (key,value) pair
if (.not. associated(keys)) then
allocate(keys)
current => keys
else
allocate(current % next)
current => current % next
end if
! Copy (key,value) pair
current % key = elem % key
current % value = elem % value
! Move to next element in bucket i
elem => elem % next
end do
end do
end function dict_keys_ii
!===============================================================================
! DICT_CLEAR Deletes and deallocates the dictionary item
!===============================================================================
subroutine dict_clear_ci(this)
class(DictCharInt) :: this
integer :: i
type(ElemKeyValueCI), pointer :: current
type(ElemKeyValueCI), pointer :: next
if (associated(this % table)) then
do i = 1, size(this % table)
current => this % table(i) % list
do while (associated(current))
if (associated(current % next)) then
next => current % next
else
nullify(next)
end if
deallocate(current)
current => next
end do
if (associated(this % table(i) % list)) &
nullify(this % table(i) % list)
end do
deallocate(this % table)
i = this % get_entry(key)
if (this % table(i) % key /= EMPTY) then
this % table(i) % key = DELETED
this % table(i) % value = EMPTY
this % entries = this % entries - 1
end if
end subroutine dict_clear_ci
end subroutine remove_ii
subroutine dict_clear_ii(this)
subroutine remove_ci(this, key)
class(DictCharInt) :: this
character(*), intent(in) :: key
integer :: i
i = this % get_entry(key)
if (this % table(i) % hash /= EMPTY) then
this % table(i) % hash = DELETED
if (allocated(this % table(i) % entry)) then
deallocate(this % table(i) % entry)
end if
this % entries = this % entries - 1
end if
end subroutine remove_ci
!===============================================================================
! NEXT_ENTRY finds the next (key,value) pair. The value of current_entry is
! updated with the (key,value) pair, and the value of i is updated with the
! index of the entry in the table. Passing in i = 0 will locate the first entry
! in the dictionary. If there are no more entries, i will be set to 0.
!===============================================================================
subroutine next_entry_ii(this, current_entry, i)
class(DictIntInt) :: this
type(DictEntryII), intent(inout) :: current_entry
integer, intent(inout) :: i
if (.not. allocated(this % table)) return
do
i = i + 1
if (i > size(this % table)) then
i = 0
exit
else if (this % table(i) % key /= EMPTY .and. &
this % table(i) % key /= DELETED) then
current_entry = this % table(i)
exit
end if
end do
end subroutine next_entry_ii
subroutine next_entry_ci(this, current_entry, i)
class(DictCharInt) :: this
type(DictEntryCI), intent(inout) :: current_entry
integer, intent(inout) :: i
if (.not. allocated(this % table)) return
do
i = i + 1
if (i > size(this % table)) then
i = 0
exit
else if (this % table(i) % hash /= EMPTY .and. &
this % table(i) % hash /= DELETED) then
current_entry = this % table(i) % entry
exit
end if
end do
end subroutine next_entry_ci
!===============================================================================
! CLEAR deletes and deallocates the dictionary item
!===============================================================================
subroutine clear_ii(this)
class(DictIntInt) :: this
if (allocated(this % table)) deallocate(this % table)
this % entries = 0
this % capacity = 0
end subroutine clear_ii
subroutine clear_ci(this)
class(DictCharInt) :: this
if (allocated(this % table)) deallocate(this % table)
this % entries = 0
this % capacity = 0
end subroutine clear_ci
!===============================================================================
! SIZE returns the number of entries in the dictionary
!===============================================================================
pure function size_ii(this) result(size)
class(DictIntInt), intent(in) :: this
integer :: size
size = this % entries
end function size_ii
pure function size_ci(this) result(size)
class(DictCharInt), intent(in) :: this
integer :: size
size = this % entries
end function size_ci
!===============================================================================
! HASH returns the hash value for a given key
!===============================================================================
pure function hash_ii(key) result(hash)
integer, intent(in) :: key
integer :: hash
hash = abs(key - 1)
end function hash_ii
pure function hash_ci(key) result(hash)
character(*), intent(in) :: key
integer :: hash
integer :: i
type(ElemKeyValueII), pointer :: current
type(ElemKeyValueII), pointer :: next
if (associated(this % table)) then
do i = 1, size(this % table)
current => this % table(i) % list
do while (associated(current))
if (associated(current % next)) then
next => current % next
else
nullify(next)
end if
deallocate(current)
current => next
end do
if (associated(this % table(i) % list)) &
nullify(this % table(i) % list)
end do
deallocate(this % table)
end if
hash = 0
end subroutine dict_clear_ii
do i = 1, len_trim(key)
hash = 31 * hash + ichar(key(i:i))
end do
hash = abs(hash - 1)
end function hash_ci
end module dict_header

12
src/geometry.F90
View file

@ -1249,8 +1249,8 @@ contains
class(Lattice), pointer :: lat ! pointer to current lattice
! Don't research places already checked
if (found(universe_dict % get_key(univ % id), map)) then
count = counts(universe_dict % get_key(univ % id), map)
if (found(universe_dict % get(univ % id), map)) then
count = counts(universe_dict % get(univ % id), map)
return
end if
@ -1258,8 +1258,8 @@ contains
! Count = 1, then quit
if (univ % id == univ_id) then
count = 1
counts(universe_dict % get_key(univ % id), map) = 1
found(universe_dict % get_key(univ % id), map) = .true.
counts(universe_dict % get(univ % id), map) = 1
found(universe_dict % get(univ % id), map) = .true.
return
end if
@ -1355,8 +1355,8 @@ contains
end if
end do
counts(universe_dict % get_key(univ % id), map) = count
found(universe_dict % get_key(univ % id), map) = .true.
counts(universe_dict % get(univ % id), map) = count
found(universe_dict % get(univ % id), map) = .true.
end function count_target

9
src/geometry_header.F90
View file

@ -374,10 +374,11 @@ contains
end if
i_material = cells(i) % material(j)
associate (mat => materials(i_material))
NUC_NAMES_LOOP: do k = 1, size(mat % names)
! Get index in nuc_temps array
i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k)))
i_nuclide = nuclide_dict % get(to_lower(mat % names(k)))
! Add temperature if it hasn't already been added
if (find(nuc_temps(i_nuclide), temperature) == -1) then
@ -388,7 +389,7 @@ contains
if (present(sab_temps) .and. mat % n_sab > 0) then
SAB_NAMES_LOOP: do k = 1, size(mat % sab_names)
! Get index in nuc_temps array
i_sab = sab_dict % get_key(to_lower(mat % sab_names(k)))
i_sab = sab_dict % get(to_lower(mat % sab_names(k)))
! Add temperature if it hasn't already been added
if (find(sab_temps(i_sab), temperature) == -1) then
@ -433,8 +434,8 @@ contains
integer(C_INT) :: err
if (allocated(cells)) then
if (cell_dict % has_key(id)) then
index = cell_dict % get_key(id)
if (cell_dict % has(id)) then
index = cell_dict % get(id)
err = 0
else
err = E_INVALID_ID

7
src/initialize.F90
View file

@ -9,7 +9,6 @@ module initialize
use bank_header, only: Bank
use constants
use dict_header, only: DictIntInt, ElemKeyValueII
use set_header, only: SetInt
use error, only: fatal_error, warning
use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,&
@ -21,7 +20,7 @@ module initialize
use message_passing
use mgxs_data, only: read_mgxs, create_macro_xs
use output, only: print_version, write_message, print_usage
use random_lcg, only: initialize_prng
use random_lcg, only: openmc_set_seed, seed
use settings
#ifdef _OPENMP
use simulation_header, only: n_threads
@ -45,6 +44,8 @@ contains
subroutine openmc_init(intracomm) bind(C)
integer, intent(in), optional :: intracomm ! MPI intracommunicator
integer :: err
! Copy the communicator to a new variable. This is done to avoid changing
! the signature of this subroutine. If MPI is being used but no communicator
! was passed, assume MPI_COMM_WORLD.
@ -83,7 +84,7 @@ contains
! Initialize random number generator -- if the user specifies a seed, it
! will be re-initialized later
call initialize_prng()
err = openmc_set_seed(seed)
! Read XML input files
call read_input_xml()

244
src/input_xml.F90
View file

@ -6,7 +6,7 @@ module input_xml
use cmfd_input, only: configure_cmfd
use cmfd_header, only: cmfd_mesh
use constants
use dict_header, only: DictIntInt, DictCharInt, ElemKeyValueCI
use dict_header, only: DictIntInt, DictCharInt, DictEntryCI
use distribution_multivariate
use distribution_univariate
use endf, only: reaction_name
@ -25,7 +25,7 @@ module input_xml
use nuclide_header
use output, only: write_message, title, header, print_plot
use plot_header
use random_lcg, only: prn, seed, initialize_prng
use random_lcg, only: prn, openmc_set_seed
use surface_header
use set_header, only: SetChar
use settings
@ -151,6 +151,7 @@ contains
integer :: temp_int
integer :: temp_int_array3(3)
integer(C_INT32_T) :: i_start, i_end
integer(C_INT64_T) :: seed
integer(C_INT) :: err
integer, allocatable :: temp_int_array(:)
integer :: n_tracks
@ -340,7 +341,7 @@ contains
! Copy random number seed if specified
if (check_for_node(root, "seed")) then
call get_node_value(root, "seed", seed)
call initialize_prng()
err = openmc_set_seed(seed)
end if
! Number of bins for logarithmic grid
@ -477,15 +478,15 @@ contains
call m % from_xml(node_mesh_list(i))
! Add mesh to dictionary
call mesh_dict % add_key(m % id, i_start + i - 1)
call mesh_dict % set(m % id, i_start + i - 1)
end associate
end do
! Shannon Entropy mesh
if (check_for_node(root, "entropy_mesh")) then
call get_node_value(root, "entropy_mesh", temp_int)
if (mesh_dict % has_key(temp_int)) then
index_entropy_mesh = mesh_dict % get_key(temp_int)
if (mesh_dict % has(temp_int)) then
index_entropy_mesh = mesh_dict % get(temp_int)
else
call fatal_error("Mesh " // to_str(temp_int) // " specified for &
&Shannon entropy does not exist.")
@ -533,8 +534,8 @@ contains
! Uniform fission source weighting mesh
if (check_for_node(root, "ufs_mesh")) then
call get_node_value(root, "ufs_mesh", temp_int)
if (mesh_dict % has_key(temp_int)) then
index_ufs_mesh = mesh_dict % get_key(temp_int)
if (mesh_dict % has(temp_int)) then
index_ufs_mesh = mesh_dict % get(temp_int)
else
call fatal_error("Mesh " // to_str(temp_int) // " specified for &
&uniform fission site method does not exist.")
@ -1044,7 +1045,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (surface_dict % has_key(s%id)) then
if (surface_dict % has(s%id)) then
call fatal_error("Two or more surfaces use the same unique ID: " &
// to_str(s%id))
end if
@ -1175,7 +1176,7 @@ contains
&"' specified on surface " // trim(to_str(s%id)))
end select
! Add surface to dictionary
call surface_dict % add_key(s%id, i)
call surface_dict % set(s%id, i)
end do
! Check to make sure a boundary condition was applied to at least one
@ -1201,7 +1202,7 @@ contains
&interior surface.")
end if
else
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
surf % i_periodic = surface_dict % get(surf % i_periodic)
end if
type is (SurfaceYPlane)
@ -1215,7 +1216,7 @@ contains
&interior surface.")
end if
else
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
surf % i_periodic = surface_dict % get(surf % i_periodic)
end if
type is (SurfaceZPlane)
@ -1229,7 +1230,7 @@ contains
&interior surface.")
end if
else
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
surf % i_periodic = surface_dict % get(surf % i_periodic)
end if
type is (SurfacePlane)
@ -1238,7 +1239,7 @@ contains
&periodic boundary condition on surface " // &
trim(to_str(surf % id)) // ".")
else
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
surf % i_periodic = surface_dict % get(surf % i_periodic)
end if
class default
@ -1315,7 +1316,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (cell_dict % has_key(c % id)) then
if (cell_dict % has(c % id)) then
call fatal_error("Two or more cells use the same unique ID: " &
// to_str(c % id))
end if
@ -1390,8 +1391,8 @@ contains
do j = 1, tokens % size()
id = tokens % data(j)
if (id < OP_UNION) then
if (surface_dict % has_key(abs(id))) then
k = surface_dict % get_key(abs(id))
if (surface_dict % has(abs(id))) then
k = surface_dict % get(abs(id))
tokens % data(j) = sign(k, id)
end if
end if
@ -1511,21 +1512,21 @@ contains
end if
! Add cell to dictionary
call cell_dict % add_key(c % id, i)
call cell_dict % set(c % id, i)
! For cells, we also need to check if there's a new universe --
! also for every cell add 1 to the count of cells for the
! specified universe
univ_id = c % universe
if (.not. cells_in_univ_dict % has_key(univ_id)) then
if (.not. cells_in_univ_dict % has(univ_id)) then
n_universes = n_universes + 1
n_cells_in_univ = 1
call universe_dict % add_key(univ_id, n_universes)
call universe_dict % set(univ_id, n_universes)
call univ_ids % push_back(univ_id)
else
n_cells_in_univ = 1 + cells_in_univ_dict % get_key(univ_id)
n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id)
end if
call cells_in_univ_dict % add_key(univ_id, n_cells_in_univ)
call cells_in_univ_dict % set(univ_id, n_cells_in_univ)
end do
@ -1559,7 +1560,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (lattice_dict % has_key(lat % id)) then
if (lattice_dict % has(lat % id)) then
call fatal_error("Two or more lattices use the same unique ID: " &
// to_str(lat % id))
end if
@ -1669,7 +1670,7 @@ contains
end if
! Add lattice to dictionary
call lattice_dict % add_key(lat % id, i)
call lattice_dict % set(lat % id, i)
end select
end do RECT_LATTICES
@ -1691,7 +1692,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (lattice_dict % has_key(lat % id)) then
if (lattice_dict % has(lat % id)) then
call fatal_error("Two or more lattices use the same unique ID: " &
// to_str(lat % id))
end if
@ -1856,7 +1857,7 @@ contains
end if
! Add lattice to dictionary
call lattice_dict % add_key(lat % id, n_rlats + i)
call lattice_dict % set(lat % id, n_rlats + i)
end select
end do HEX_LATTICES
@ -1871,7 +1872,7 @@ contains
u % id = univ_ids % data(i)
! Allocate cell list
n_cells_in_univ = cells_in_univ_dict % get_key(u % id)
n_cells_in_univ = cells_in_univ_dict % get(u % id)
allocate(u % cells(n_cells_in_univ))
u % cells(:) = 0
@ -1890,7 +1891,7 @@ contains
do i = 1, n_cells
! Get index in universes array
j = universe_dict % get_key(cells(i) % universe)
j = universe_dict % get(cells(i) % universe)
! Set the first zero entry in the universe cells array to the index in the
! global cells array
@ -2007,14 +2008,14 @@ contains
! Creating dictionary that maps the name of the material to the entry
do i = 1, size(libraries)
do j = 1, size(libraries(i) % materials)
call library_dict % add_key(to_lower(libraries(i) % materials(j)), i)
call library_dict % set(to_lower(libraries(i) % materials(j)), i)
end do
end do
! Check that 0K nuclides are listed in the cross_sections.xml file
if (allocated(res_scat_nuclides)) then
do i = 1, size(res_scat_nuclides)
if (.not. library_dict % has_key(to_lower(res_scat_nuclides(i)))) then
if (.not. library_dict % has(to_lower(res_scat_nuclides(i)))) then
call fatal_error("Could not find resonant scatterer " &
// trim(res_scat_nuclides(i)) // " in cross_sections.xml file!")
end if
@ -2102,7 +2103,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (material_dict % has_key(mat % id)) then
if (material_dict % has(mat % id)) then
call fatal_error("Two or more materials use the same unique ID: " &
// to_str(mat % id))
end if
@ -2311,11 +2312,11 @@ contains
ALL_NUCLIDES: do j = 1, mat % n_nuclides
! Check that this nuclide is listed in the cross_sections.xml file
name = trim(names % data(j))
if (.not. library_dict % has_key(to_lower(name))) then
if (.not. library_dict % has(to_lower(name))) then
call fatal_error("Could not find nuclide " // trim(name) &
// " in cross_sections data file!")
end if
i_library = library_dict % get_key(to_lower(name))
i_library = library_dict % get(to_lower(name))
if (run_CE) then
! Check to make sure cross-section is continuous energy neutron table
@ -2327,13 +2328,13 @@ contains
! If this nuclide hasn't been encountered yet, we need to add its name
! and alias to the nuclide_dict
if (.not. nuclide_dict % has_key(to_lower(name))) then
if (.not. nuclide_dict % has(to_lower(name))) then
index_nuclide = index_nuclide + 1
mat % nuclide(j) = index_nuclide
call nuclide_dict % add_key(to_lower(name), index_nuclide)
call nuclide_dict % set(to_lower(name), index_nuclide)
else
mat % nuclide(j) = nuclide_dict % get_key(to_lower(name))
mat % nuclide(j) = nuclide_dict % get(to_lower(name))
end if
! Copy name and atom/weight percent
@ -2406,14 +2407,14 @@ contains
end if
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. library_dict % has_key(to_lower(name))) then
if (.not. library_dict % has(to_lower(name))) then
call fatal_error("Could not find S(a,b) table " // trim(name) &
// " in cross_sections.xml file!")
end if
! Find index in xs_listing and set the name and alias according to the
! listing
i_library = library_dict % get_key(to_lower(name))
i_library = library_dict % get(to_lower(name))
if (run_CE) then
! Check to make sure cross-section is continuous energy neutron table
@ -2425,19 +2426,19 @@ contains
! If this S(a,b) table hasn't been encountered yet, we need to add its
! name and alias to the sab_dict
if (.not. sab_dict % has_key(to_lower(name))) then
if (.not. sab_dict % has(to_lower(name))) then
index_sab = index_sab + 1
mat % i_sab_tables(j) = index_sab
call sab_dict % add_key(to_lower(name), index_sab)
call sab_dict % set(to_lower(name), index_sab)
else
mat % i_sab_tables(j) = sab_dict % get_key(to_lower(name))
mat % i_sab_tables(j) = sab_dict % get(to_lower(name))
end if
end do
end if
end if
! Add material to dictionary
call material_dict % add_key(mat % id, i)
call material_dict % set(mat % id, i)
end do
! Set total number of nuclides and S(a,b) tables
@ -2463,6 +2464,7 @@ contains
integer :: filter_id ! user-specified identifier for filter
integer :: i_filt ! index in filters array
integer :: i_filter_mesh ! index of mesh filter
integer :: i_elem ! index of entry in dictionary
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
integer :: n_filter ! number of filters
@ -2500,8 +2502,7 @@ contains
type(XMLNode), allocatable :: node_filt_list(:)
type(XMLNode), allocatable :: node_trigger_list(:)
type(XMLNode), allocatable :: node_deriv_list(:)
type(ElemKeyValueCI), pointer :: scores
type(ElemKeyValueCI), pointer :: next
type(DictEntryCI) :: elem
! Check if tallies.xml exists
filename = trim(path_input) // "tallies.xml"
@ -2557,7 +2558,7 @@ contains
call m % from_xml(node_mesh_list(i))
! Add mesh to dictionary
call mesh_dict % add_key(m % id, i_start + i - 1)
call mesh_dict % set(m % id, i_start + i - 1)
end do
! We only need the mesh info for plotting
@ -2586,7 +2587,7 @@ contains
call tally_derivs(i) % from_xml(node_deriv_list(i))
! Update tally derivative dictionary
call tally_deriv_dict % add_key(tally_derivs(i) % id, i)
call tally_deriv_dict % set(tally_derivs(i) % id, i)
end do
! ==========================================================================
@ -2612,7 +2613,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (filter_dict % has_key(filter_id)) then
if (filter_dict % has(filter_id)) then
call fatal_error("Two or more filters use the same unique ID: " &
// to_str(filter_id))
end if
@ -2682,7 +2683,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (tally_dict % has_key(t % id)) then
if (tally_dict % has(t % id)) then
call fatal_error("Two or more tallies use the same unique ID: " &
// to_str(t % id))
end if
@ -2717,8 +2718,8 @@ contains
do j = 1, n_filter
! Get pointer to filter
if (filter_dict % has_key(temp_filter(j))) then
i_filt = filter_dict % get_key(temp_filter(j))
if (filter_dict % has(temp_filter(j))) then
i_filt = filter_dict % get(temp_filter(j))
f => filters(i_filt)
else
call fatal_error("Could not find filter " &
@ -2774,14 +2775,14 @@ contains
word = to_lower(sarray(j))
! Search through nuclides
if (.not. nuclide_dict % has_key(word)) then
if (.not. nuclide_dict % has(word)) then
call fatal_error("Could not find the nuclide " &
// trim(word) // " specified in tally " &
// trim(to_str(t % id)) // " in any material.")
end if
! Set bin to index in nuclides array
t % nuclide_bins(j) = nuclide_dict % get_key(word)
t % nuclide_bins(j) = nuclide_dict % get(word)
end do
! Set number of nuclide bins
@ -2819,7 +2820,7 @@ contains
score_name = trim(sarray(j))
! Append the score to the list of possible trigger scores
if (trigger_on) call trigger_scores % add_key(trim(score_name), j)
if (trigger_on) call trigger_scores % set(trim(score_name), j)
do imomstr = 1, size(MOMENT_STRS)
if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
@ -3145,7 +3146,9 @@ contains
filters(i_filt_start) % obj % n_bins = product(m % dimension + 1)
! Set ID
err = openmc_filter_set_id(i_filt_start, i_filt_start)
call openmc_get_filter_next_id(filter_id)
err = openmc_filter_set_id(i_filt_start, filter_id)
! Add surface filter
allocate(SurfaceFilter :: filters(i_filt_end) % obj)
@ -3166,7 +3169,8 @@ contains
filt % current = .true.
! Set ID
err = openmc_filter_set_id(i_filt_end, i_filt_end)
call openmc_get_filter_next_id(filter_id)
err = openmc_filter_set_id(i_filt_end, filter_id)
end select
! Copy filter indices to resized array
@ -3399,15 +3403,7 @@ contains
score_name = trim(to_lower(sarray(j)))
if (score_name == "all") then
scores => trigger_scores % keys()
do while (associated(scores))
next => scores % next
deallocate(scores)
scores => next
t % n_triggers = t % n_triggers + 1
end do
t % n_triggers = t % n_triggers + trigger_scores % size()
else
t % n_triggers = t % n_triggers + 1
end if
@ -3466,16 +3462,19 @@ contains
! Expand "all" to include TriggerObjects for each score in tally
if (score_name == "all") then
scores => trigger_scores % keys()
! Loop over all tally scores
do while (associated(scores))
score_name = trim(scores % key)
i_elem = 0
do
! Move to next score
call trigger_scores % next_entry(elem, i_elem)
if (i_elem == 0) exit
score_name = trim(elem % key)
! Store the score name and index in the trigger
t % triggers(trig_ind) % score_name = trim(score_name)
t % triggers(trig_ind) % score_index = &
trigger_scores % get_key(trim(score_name))
t % triggers(trig_ind) % score_name = score_name
t % triggers(trig_ind) % score_index = elem % value
! Set the trigger convergence threshold type
select case (temp_str)
@ -3493,11 +3492,6 @@ contains
! Store the trigger convergence threshold
t % triggers(trig_ind) % threshold = threshold
! Move to next score
next => scores % next
deallocate(scores)
scores => next
! Increment the overall trigger index
trig_ind = trig_ind + 1
end do
@ -3508,7 +3502,7 @@ contains
! Store the score name and index
t % triggers(trig_ind) % score_name = trim(score_name)
t % triggers(trig_ind) % score_index = &
trigger_scores % get_key(trim(score_name))
trigger_scores % get(trim(score_name))
! Check if an invalid score was set for the trigger
if (t % triggers(trig_ind) % score_index == 0) then
@ -3585,7 +3579,7 @@ contains
end if
! Add tally to dictionary
call tally_dict % add_key(t % id, i)
call tally_dict % set(t % id, i)
end associate
end do READ_TALLIES
@ -3657,7 +3651,7 @@ contains
end if
! Check to make sure 'id' hasn't been used
if (plot_dict % has_key(pl % id)) then
if (plot_dict % has(pl % id)) then
call fatal_error("Two or more plots use the same unique ID: " &
// to_str(pl % id))
end if
@ -3842,8 +3836,8 @@ contains
! Add RGB
if (pl % color_by == PLOT_COLOR_CELLS) then
if (cell_dict % has_key(col_id)) then
col_id = cell_dict % get_key(col_id)
if (cell_dict % has(col_id)) then
col_id = cell_dict % get(col_id)
call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
else
call fatal_error("Could not find cell " // trim(to_str(col_id)) &
@ -3852,8 +3846,8 @@ contains
else if (pl % color_by == PLOT_COLOR_MATS) then
if (material_dict % has_key(col_id)) then
col_id = material_dict % get_key(col_id)
if (material_dict % has(col_id)) then
col_id = material_dict % get(col_id)
call get_node_array(node_col, "rgb", pl % colors(col_id) % rgb)
else
call fatal_error("Could not find material " &
@ -3957,8 +3951,8 @@ contains
end if
! Check if the specified tally mesh exists
if (mesh_dict % has_key(meshid)) then
pl % meshlines_mesh => meshes(mesh_dict % get_key(meshid))
if (mesh_dict % has(meshid)) then
pl % meshlines_mesh => meshes(mesh_dict % get(meshid))
if (meshes(meshid) % type /= LATTICE_RECT) then
call fatal_error("Non-rectangular mesh specified in &
&meshlines for plot " // trim(to_str(pl % id)))
@ -4017,8 +4011,8 @@ contains
if (pl % color_by == PLOT_COLOR_CELLS) then
if (cell_dict % has_key(col_id)) then
iarray(j) = cell_dict % get_key(col_id)
if (cell_dict % has(col_id)) then
iarray(j) = cell_dict % get(col_id)
else
call fatal_error("Could not find cell " &
// trim(to_str(col_id)) // " specified in the mask in &
@ -4027,8 +4021,8 @@ contains
else if (pl % color_by == PLOT_COLOR_MATS) then
if (material_dict % has_key(col_id)) then
iarray(j) = material_dict % get_key(col_id)
if (material_dict % has(col_id)) then
iarray(j) = material_dict % get(col_id)
else
call fatal_error("Could not find material " &
// trim(to_str(col_id)) // " specified in the mask in &
@ -4056,7 +4050,7 @@ contains
end if
! Add plot to dictionary
call plot_dict % add_key(pl % id, i)
call plot_dict % set(pl % id, i)
end do READ_PLOTS
@ -4437,8 +4431,8 @@ contains
name = materials(i) % names(j)
if (.not. already_read % contains(name)) then
i_library = library_dict % get_key(to_lower(name))
i_nuclide = nuclide_dict % get_key(to_lower(name))
i_library = library_dict % get(to_lower(name))
i_nuclide = nuclide_dict % get(to_lower(name))
call write_message('Reading ' // trim(name) // ' from ' // &
trim(libraries(i_library) % path), 6)
@ -4497,8 +4491,8 @@ contains
name = materials(i) % sab_names(j)
if (.not. already_read % contains(name)) then
i_library = library_dict % get_key(to_lower(name))
i_sab = sab_dict % get_key(to_lower(name))
i_library = library_dict % get(to_lower(name))
i_sab = sab_dict % get(to_lower(name))
call write_message('Reading ' // trim(name) // ' from ' // &
trim(libraries(i_library) % path), 6)
@ -4663,12 +4657,12 @@ contains
! ADJUST UNIVERSE INDEX FOR EACH CELL
associate (c => cells(i))
id = c%universe
if (universe_dict%has_key(id)) then
c%universe = universe_dict%get_key(id)
id = c % universe
if (universe_dict % has(id)) then
c % universe = universe_dict % get(id)
else
call fatal_error("Could not find universe " // trim(to_str(id)) &
&// " specified on cell " // trim(to_str(c%id)))
&// " specified on cell " // trim(to_str(c % id)))
end if
! =======================================================================
@ -4676,11 +4670,11 @@ contains
if (c % material(1) == NONE) then
id = c % fill
if (universe_dict % has_key(id)) then
if (universe_dict % has(id)) then
c % type = FILL_UNIVERSE
c % fill = universe_dict % get_key(id)
elseif (lattice_dict % has_key(id)) then
lid = lattice_dict % get_key(id)
c % fill = universe_dict % get(id)
elseif (lattice_dict % has(id)) then
lid = lattice_dict % get(id)
c % type = FILL_LATTICE
c % fill = lid
else
@ -4693,9 +4687,9 @@ contains
id = c % material(j)
if (id == MATERIAL_VOID) then
c % type = FILL_MATERIAL
else if (material_dict % has_key(id)) then
else if (material_dict % has(id)) then
c % type = FILL_MATERIAL
c % material(j) = material_dict % get_key(id)
c % material(j) = material_dict % get(id)
else
call fatal_error("Could not find material " // trim(to_str(id)) &
// " specified on cell " // trim(to_str(c % id)))
@ -4709,41 +4703,41 @@ contains
! ADJUST UNIVERSE INDICES FOR EACH LATTICE
do i = 1, n_lattices
lat => lattices(i)%obj
lat => lattices(i) % obj
select type (lat)
type is (RectLattice)
do m = 1, lat%n_cells(3)
do k = 1, lat%n_cells(2)
do j = 1, lat%n_cells(1)
id = lat%universes(j,k,m)
if (universe_dict%has_key(id)) then
lat%universes(j,k,m) = universe_dict%get_key(id)
do m = 1, lat % n_cells(3)
do k = 1, lat % n_cells(2)
do j = 1, lat % n_cells(1)
id = lat % universes(j,k,m)
if (universe_dict % has(id)) then
lat % universes(j,k,m) = universe_dict % get(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat%id)))
&// trim(to_str(lat % id)))
end if
end do
end do
end do
type is (HexLattice)
do m = 1, lat%n_axial
do k = 1, 2*lat%n_rings - 1
do j = 1, 2*lat%n_rings - 1
if (j + k < lat%n_rings + 1) then
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
if (j + k < lat % n_rings + 1) then
cycle
else if (j + k > 3*lat%n_rings - 1) then
else if (j + k > 3*lat % n_rings - 1) then
cycle
end if
id = lat%universes(j, k, m)
if (universe_dict%has_key(id)) then
lat%universes(j, k, m) = universe_dict%get_key(id)
id = lat % universes(j, k, m)
if (universe_dict % has(id)) then
lat % universes(j, k, m) = universe_dict % get(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat%id)))
&// trim(to_str(lat % id)))
end if
end do
end do
@ -4751,13 +4745,13 @@ contains
end select
if (lat%outer /= NO_OUTER_UNIVERSE) then
if (universe_dict%has_key(lat%outer)) then
lat%outer = universe_dict%get_key(lat%outer)
if (lat % outer /= NO_OUTER_UNIVERSE) then
if (universe_dict % has(lat % outer)) then
lat % outer = universe_dict % get(lat % outer)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(lat%outer)) &
&// " specified on lattice " // trim(to_str(lat%id)))
&// trim(to_str(lat % outer)) &
&// " specified on lattice " // trim(to_str(lat % id)))
end if
end if

12
src/material_header.F90
View file

@ -277,8 +277,8 @@ contains
integer(C_INT) :: err
if (allocated(materials)) then
if (material_dict % has_key(id)) then
index = material_dict % get_key(id)
if (material_dict % has(id)) then
index = material_dict % get(id)
err = 0
else
err = E_INVALID_ID
@ -345,7 +345,7 @@ contains
call move_alloc(FROM=new_p0, TO=m % p0)
! Append new nuclide/density
k = nuclide_dict % get_key(to_lower(name_))
k = nuclide_dict % get(to_lower(name_))
m % nuclide(n + 1) = k
m % atom_density(n + 1) = density
m % density = m % density + density
@ -415,7 +415,7 @@ contains
if (index >= 1 .and. index <= n_materials) then
materials(index) % id = id
call material_dict % add_key(id, index)
call material_dict % set(id, index)
err = 0
else
err = E_OUT_OF_BOUNDS
@ -469,12 +469,12 @@ contains
call c_f_pointer(name(i), string, [10])
name_ = to_lower(to_f_string(string))
if (.not. nuclide_dict % has_key(name_)) then
if (.not. nuclide_dict % has(name_)) then
err = openmc_load_nuclide(string)
if (err < 0) return
end if
m % nuclide(i) = nuclide_dict % get_key(name_)
m % nuclide(i) = nuclide_dict % get(name_)
m % atom_density(i) = density(i)
end do
m % n_nuclides = n

2
src/mesh_header.F90
View file

@ -61,7 +61,7 @@ contains
call get_node_value(node, "id", this % id)
! Check to make sure 'id' hasn't been used
if (mesh_dict % has_key(this % id)) then
if (mesh_dict % has(this % id)) then
call fatal_error("Two or more meshes use the same unique ID: " &
// to_str(this % id))
end if

3
src/mgxs_data.F90
View file

@ -26,7 +26,6 @@ contains
subroutine read_mgxs()
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: i_xsdata ! index in <xsdata> list
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
integer :: representation ! Data representation
@ -36,7 +35,6 @@ contains
integer(HID_T) :: file_id
integer(HID_T) :: xsdata_group
logical :: file_exists
type(DictCharInt) :: xsdata_dict
type(VectorReal), allocatable :: temps(:)
character(MAX_WORD_LEN) :: word
integer, allocatable :: array(:)
@ -86,7 +84,6 @@ contains
name = mat % names(j)
if (.not. already_read % contains(name)) then
i_xsdata = xsdata_dict % get_key(to_lower(name))
i_nuclide = mat % nuclide(j)
call write_message("Loading " // trim(name) // " data...", 6)

14
src/nuclide_header.F90
View file

@ -581,7 +581,7 @@ contains
do i = 1, size(this % reactions)
call MTs % push_back(this % reactions(i) % MT)
call this % reaction_index % add_key(this % reactions(i) % MT, i)
call this % reaction_index % set(this % reactions(i) % MT, i)
associate (rx => this % reactions(i))
! Skip total inelastic level scattering, gas production cross sections
@ -862,8 +862,8 @@ contains
name_ = to_f_string(name)
if (allocated(nuclides)) then
if (nuclide_dict % has_key(to_lower(name_))) then
index = nuclide_dict % get_key(to_lower(name_))
if (nuclide_dict % has(to_lower(name_))) then
index = nuclide_dict % get(to_lower(name_))
err = 0
else
err = E_DATA
@ -895,8 +895,8 @@ contains
name_ = to_f_string(name)
err = 0
if (.not. nuclide_dict % has_key(to_lower(name_))) then
if (library_dict % has_key(to_lower(name_))) then
if (.not. nuclide_dict % has(to_lower(name_))) then
if (library_dict % has(to_lower(name_))) then
! allocate extra space in nuclides array
n = n_nuclides
allocate(new_nuclides(n + 1))
@ -904,7 +904,7 @@ contains
call move_alloc(FROM=new_nuclides, TO=nuclides)
n = n + 1
i_library = library_dict % get_key(to_lower(name_))
i_library = library_dict % get(to_lower(name_))
! Open file and make sure version is sufficient
file_id = file_open(libraries(i_library) % path, 'r')
@ -919,7 +919,7 @@ contains
call file_close(file_id)
! Add entry to nuclide dictionary
call nuclide_dict % add_key(to_lower(name_), n)
call nuclide_dict % set(to_lower(name_), n)
n_nuclides = n
! Assign resonant scattering data

109
src/random_lcg.F90
View file

@ -1,5 +1,7 @@
module random_lcg
use, intrinsic :: ISO_C_BINDING
use constants
implicit none
@ -7,28 +9,29 @@ module random_lcg
private
save
! Random number seed
integer(8), public :: seed = 1_8
! Starting seed
integer(C_INT64_T), public, bind(C) :: seed = 1_8
integer(8) :: prn_seed0 ! original seed
integer(8) :: prn_seed(N_STREAMS) ! current seed
integer(8) :: prn_mult ! multiplication factor, g
integer(8) :: prn_add ! additive factor, c
integer :: prn_bits ! number of bits, M
integer(8) :: prn_mod ! 2^M
integer(8) :: prn_mask ! 2^M - 1
integer(8) :: prn_stride ! stride between particles
real(8) :: prn_norm ! 2^(-M)
integer :: stream ! current RNG stream
! LCG parameters
integer(C_INT64_T), parameter :: prn_mult = 2806196910506780709_8 ! multiplication factor, g
integer(C_INT64_T), parameter :: prn_add = 1_8 ! additive factor, c
integer, parameter :: prn_bits = 63 ! number of bits, M
integer(C_INT64_T), parameter :: prn_mod = ibset(0_8, prn_bits) ! 2^M
integer(C_INT64_T), parameter :: prn_mask = not(prn_mod) ! 2^M - 1
integer(C_INT64_T), parameter :: prn_stride = 152917_8 ! stride between particles
real(C_DOUBLE), parameter :: prn_norm = 2._8**(-prn_bits) ! 2^(-M)
! Current PRNG state
integer(C_INT64_T) :: prn_seed(N_STREAMS) ! current seed
integer :: stream ! current RNG stream
!$omp threadprivate(prn_seed, stream)
public :: prn
public :: future_prn
public :: initialize_prng
public :: set_particle_seed
public :: advance_prn_seed
public :: prn_set_stream
public :: openmc_set_seed
contains
@ -38,10 +41,10 @@ contains
function prn() result(pseudo_rn)
real(8) :: pseudo_rn
real(C_DOUBLE) :: pseudo_rn
! This algorithm uses bit-masking to find the next integer(8) value to be
! used to calculate the random number
! This algorithm uses bit-masking to find the next integer(C_INT64_T) value
! to be used to calculate the random number
prn_seed(stream) = iand(prn_mult*prn_seed(stream) + prn_add, prn_mask)
@ -59,40 +62,14 @@ contains
function future_prn(n) result(pseudo_rn)
integer(8), intent(in) :: n ! number of prns to skip
integer(C_INT64_T), intent(in) :: n ! number of prns to skip
real(8) :: pseudo_rn
real(C_DOUBLE) :: pseudo_rn
pseudo_rn = future_seed(n, prn_seed(stream)) * prn_norm
end function future_prn
!===============================================================================
! INITIALIZE_PRNG sets up the random number generator, determining the seed and
! values for g, c, and m.
!===============================================================================
subroutine initialize_prng()
integer :: i
prn_seed0 = seed
!$omp parallel
do i = 1, N_STREAMS
prn_seed(i) = prn_seed0 + i - 1
end do
stream = STREAM_TRACKING
!$omp end parallel
prn_mult = 2806196910506780709_8
prn_add = 1_8
prn_bits = 63
prn_mod = ibset(0_8, prn_bits) ! clever way of calculating 2**bits
prn_mask = prn_mod - 1_8
prn_stride = 152917_8
prn_norm = 2._8**(-prn_bits)
end subroutine initialize_prng
!===============================================================================
! SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the
! particle
@ -100,12 +77,12 @@ contains
subroutine set_particle_seed(id)
integer(8), intent(in) :: id
integer(C_INT64_T), intent(in) :: id
integer :: i
do i = 1, N_STREAMS
prn_seed(i) = future_seed(id*prn_stride, prn_seed0 + i - 1)
prn_seed(i) = future_seed(id*prn_stride, seed + i - 1)
end do
end subroutine set_particle_seed
@ -117,7 +94,7 @@ contains
subroutine advance_prn_seed(n)
integer(8), intent(in) :: n ! number of seeds to skip
integer(C_INT64_T), intent(in) :: n ! number of seeds to skip
prn_seed(stream) = future_seed(n, prn_seed(stream))
@ -132,15 +109,15 @@ contains
function future_seed(n, seed) result(new_seed)
integer(8), intent(in) :: n ! number of seeds to skip
integer(8), intent(in) :: seed ! original seed
integer(8) :: new_seed ! new seed
integer(C_INT64_T), intent(in) :: n ! number of seeds to skip
integer(C_INT64_T), intent(in) :: seed ! original seed
integer(C_INT64_T) :: new_seed ! new seed
integer(8) :: nskip ! positive number of seeds to skip
integer(8) :: g ! original multiplicative constant
integer(8) :: c ! original additive constnat
integer(8) :: g_new ! new effective multiplicative constant
integer(8) :: c_new ! new effective additive constant
integer(C_INT64_T) :: nskip ! positive number of seeds to skip
integer(C_INT64_T) :: g ! original multiplicative constant
integer(C_INT64_T) :: c ! original additive constnat
integer(C_INT64_T) :: g_new ! new effective multiplicative constant
integer(C_INT64_T) :: c_new ! new effective additive constant
! In cases where we want to skip backwards, we add the period of the random
! number generator until the number of PRNs to skip is positive since
@ -198,4 +175,26 @@ contains
end subroutine prn_set_stream
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_set_seed(new_seed) result(err) bind(C)
! Saves the starting seed and sets up the PRNG thread state
integer(C_INT64_T), value, intent(in) :: new_seed
integer(C_INT) :: err
integer :: i
err = 0
seed = new_seed
!$omp parallel
do i = 1, N_STREAMS
prn_seed(i) = seed + i - 1
end do
stream = STREAM_TRACKING
!$omp end parallel
end function openmc_set_seed
end module random_lcg

20
src/tallies/tally.F90
View file

@ -5,6 +5,7 @@ module tally
use algorithm, only: binary_search
use constants
use cross_section, only: multipole_deriv_eval
use dict_header, only: EMPTY
use error, only: fatal_error
use geometry_header
use math, only: t_percentile, calc_pn, calc_rn
@ -246,8 +247,7 @@ contains
! multiplicities of one.
score = p % last_wgt * flux
else
m = nuclides(p % event_nuclide) % reaction_index % &
get_key(p % event_MT)
m = nuclides(p % event_nuclide) % reaction_index % get(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p % event_nuclide) % reactions(m))
@ -273,8 +273,7 @@ contains
! multiplicities of one.
score = p % last_wgt * flux
else
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
m = nuclides(p % event_nuclide) % reaction_index % get(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p % event_nuclide) % reactions(m))
@ -300,8 +299,7 @@ contains
! multiplicities of one.
score = p % last_wgt * flux
else
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
m = nuclides(p % event_nuclide) % reaction_index % get(p % event_MT)
! Get yield and apply to score
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
@ -1150,9 +1148,8 @@ contains
score = ZERO
if (i_nuclide > 0) then
if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then
m = nuclides(i_nuclide)%reaction_index%get_key(score_bin)
m = nuclides(i_nuclide) % reaction_index % get(score_bin)
if (m /= EMPTY) then
! Retrieve temperature and energy grid index and interpolation
! factor
i_temp = micro_xs(i_nuclide) % index_temp
@ -1190,9 +1187,8 @@ contains
! Get index in nuclides array
i_nuc = materials(p % material) % nuclide(l)
if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then
m = nuclides(i_nuc)%reaction_index%get_key(score_bin)
m = nuclides(i_nuc) % reaction_index % get(score_bin)
if (m /= EMPTY) then
! Retrieve temperature and energy grid index and
! interpolation factor
i_temp = micro_xs(i_nuc) % index_temp

10
src/tallies/tally_derivative_header.F90
View file

@ -1,7 +1,7 @@
module tally_derivative_header
use constants
use dict_header, only: DictIntInt
use dict_header, only: DictIntInt, EMPTY
use error, only: fatal_error
use nuclide_header, only: nuclide_dict
use string, only: to_str, to_lower
@ -40,6 +40,7 @@ contains
character(MAX_WORD_LEN) :: temp_str
character(MAX_WORD_LEN) :: word
integer :: val
! Copy the derivative id.
if (check_for_node(node, "id")) then
@ -56,7 +57,7 @@ contains
end if
! Make sure this id has not already been used.
if (tally_deriv_dict % has_key(this % id)) then
if (tally_deriv_dict % has(this % id)) then
call fatal_error("Two or more <derivative>'s use the same unique &
&ID: " // trim(to_str(this % id)))
end if
@ -74,12 +75,13 @@ contains
call get_node_value(node, "nuclide", word)
word = trim(to_lower(word))
if (.not. nuclide_dict % has_key(word)) then
val = nuclide_dict % get(word)
if (val == EMPTY) then
call fatal_error("Could not find the nuclide " &
// trim(word) // " specified in derivative " &
// trim(to_str(this % id)) // " in any material.")
end if
this % diff_nuclide = nuclide_dict % get_key(word)
this % diff_nuclide = val
case("temperature")
this % variable = DIFF_TEMPERATURE

15
src/tallies/tally_filter_cell.F90
View file

@ -5,6 +5,7 @@ module tally_filter_cell
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error
use hdf5_interface
use geometry_header
@ -55,11 +56,13 @@ contains
type(TallyFilterMatch), intent(inout) :: match
integer :: i
integer :: val
! Iterate over coordinate levels to see with cells match
do i = 1, p % n_coord
if (this % map % has_key(p % coord(i) % cell)) then
call match % bins % push_back(this % map % get_key(p % coord(i) % cell))
val = this % map % get(p % coord(i) % cell)
if (val /= EMPTY) then
call match % bins % push_back(val)
call match % weights % push_back(ONE)
end if
end do
@ -87,12 +90,14 @@ contains
class(CellFilter), intent(inout) :: this
integer :: i, id
integer :: val
! Convert ids to indices.
do i = 1, this % n_bins
id = this % cells(i)
if (cell_dict % has_key(id)) then
this % cells(i) = cell_dict % get_key(id)
val = cell_dict % get(id)
if (val /= EMPTY) then
this % cells(i) = val
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
&// " specified on tally filter.")
@ -101,7 +106,7 @@ contains
! Generate mapping from cell indices to filter bins.
do i = 1, this % n_bins
call this % map % add_key(this % cells(i), i)
call this % map % set(this % cells(i), i)
end do
end subroutine initialize_cell

16
src/tallies/tally_filter_cellborn.F90
View file

@ -5,6 +5,7 @@ module tally_filter_cellborn
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error
use hdf5_interface
use geometry_header
@ -54,8 +55,11 @@ contains
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
if (this % map % has_key(p % cell_born)) then
call match % bins % push_back(this % map % get_key(p % cell_born))
integer :: val
val = this % map % get(p % cell_born)
if (val /= EMPTY) then
call match % bins % push_back(val)
call match % weights % push_back(ONE)
end if
@ -81,12 +85,14 @@ contains
class(CellbornFilter), intent(inout) :: this
integer :: i, id
integer :: val
! Convert ids to indices.
do i = 1, this % n_bins
id = this % cells(i)
if (cell_dict % has_key(id)) then
this % cells(i) = cell_dict % get_key(id)
val = cell_dict % get(id)
if (val /= EMPTY) then
this % cells(i) = val
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
&// " specified on tally filter.")
@ -95,7 +101,7 @@ contains
! Generate mapping from cell indices to filter bins.
do i = 1, this % n_bins
call this % map % add_key(this % cells(i), i)
call this % map % set(this % cells(i), i)
end do
end subroutine initialize_cellborn

7
src/tallies/tally_filter_cellfrom.F90
View file

@ -5,6 +5,7 @@ module tally_filter_cellfrom
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error
use hdf5_interface
use geometry_header
@ -39,11 +40,13 @@ contains
type(TallyFilterMatch), intent(inout) :: match
integer :: i
integer :: val
! Starting one coordinate level deeper, find the next bin.
do i = 1, p % last_n_coord
if (this % map % has_key(p % last_cell(i))) then
call match % bins % push_back(this % map % get_key(p % last_cell(i)))
val = this % map % get(p % last_cell(i))
if (val /= EMPTY) then
call match % bins % push_back(val)
call match % weights % push_back(ONE)
exit
end if

9
src/tallies/tally_filter_distribcell.F90
View file

@ -5,6 +5,7 @@ module tally_filter_distribcell
use hdf5, only: HID_T
use constants
use dict_header, only: EMPTY
use error
use geometry_header
use hdf5_interface
@ -96,11 +97,13 @@ contains
class(DistribcellFilter), intent(inout) :: this
integer :: id
integer :: val
! Convert id to index.
id = this % cell
if (cell_dict % has_key(id)) then
this % cell = cell_dict % get_key(id)
val = cell_dict % get(id)
if (val /= EMPTY) then
this % cell = val
this % n_bins = cells(this % cell) % instances
else
call fatal_error("Could not find cell " // trim(to_str(id)) &
@ -159,7 +162,7 @@ contains
n = size(univ % cells)
! Write to the geometry stack
i_univ = universe_dict % get_key(univ % id)
i_univ = universe_dict % get(univ % id)
if (i_univ == root_universe) then
path = trim(path) // "u" // to_str(univ%id)
else

21
src/tallies/tally_filter_header.F90
View file

@ -19,6 +19,7 @@ module tally_filter_header
public :: openmc_filter_get_id
public :: openmc_filter_set_id
public :: openmc_get_filter_index
public :: openmc_get_filter_next_id
!===============================================================================
! TALLYFILTERMATCH stores every valid bin and weight for a filter
@ -122,6 +123,10 @@ module tally_filter_header
! Dictionary that maps user IDs to indices in 'filters'
type(DictIntInt), public :: filter_dict
! The largest filter ID that has been specified in the system. This is useful
! in case the code needs to find an ID for a new filter.
integer :: largest_filter_id
contains
!===============================================================================
@ -140,6 +145,7 @@ contains
n_filters = 0
if (allocated(filters)) deallocate(filters)
call filter_dict % clear()
largest_filter_id = 0
end subroutine free_memory_tally_filter
!===============================================================================
@ -204,7 +210,8 @@ contains
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
filters(index) % obj % id = id
call filter_dict % add_key(id, index)
call filter_dict % set(id, index)
if (id > largest_filter_id) largest_filter_id = id
err = 0
else
@ -225,8 +232,8 @@ contains
integer(C_INT) :: err
if (allocated(filters)) then
if (filter_dict % has_key(id)) then
index = filter_dict % get_key(id)
if (filter_dict % has(id)) then
index = filter_dict % get(id)
err = 0
else
err = E_INVALID_ID
@ -238,4 +245,12 @@ contains
end if
end function openmc_get_filter_index
subroutine openmc_get_filter_next_id(id) bind(C)
! Returns an ID number that has not been used by any other filters.
integer(C_INT32_T), intent(out) :: id
id = largest_filter_id + 1
end subroutine openmc_get_filter_next_id
end module tally_filter_header

23
src/tallies/tally_filter_material.F90
View file

@ -5,7 +5,7 @@ module tally_filter_material
use hdf5, only: HID_T
use constants
use dict_header, only: DictIntInt
use dict_header, only: DictIntInt, EMPTY
use error
use hdf5_interface
use material_header, only: materials, material_dict
@ -56,10 +56,13 @@ contains
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
if (this % map % has_key(p % material)) then
call match % bins % push_back(this % map % get_key(p % material))
call match % weights % push_back(ONE)
end if
integer :: val
val = this % map % get(p % material)
if (val /= EMPTY) then
call match % bins % push_back(val)
call match % weights % push_back(ONE)
end if
end subroutine get_all_bins_material
@ -84,12 +87,14 @@ contains
class(MaterialFilter), intent(inout) :: this
integer :: i, id
integer :: val
! Convert ids to indices.
do i = 1, this % n_bins
id = this % materials(i)
if (material_dict % has_key(id)) then
this % materials(i) = material_dict % get_key(id)
val = material_dict % get(id)
if (val /= EMPTY) then
this % materials(i) = val
else
call fatal_error("Could not find material " // trim(to_str(id)) &
&// " specified on a tally filter.")
@ -98,7 +103,7 @@ contains
! Generate mapping from material indices to filter bins.
do i = 1, this % n_bins
call this % map % add_key(this % materials(i), i)
call this % map % set(this % materials(i), i)
end do
end subroutine initialize_material
@ -166,7 +171,7 @@ contains
! Generate mapping from material indices to filter bins.
call f % map % clear()
do i = 1, n
call f % map % add_key(f % materials(i), i)
call f % map % set(f % materials(i), i)
end do
class default

7
src/tallies/tally_filter_mesh.F90
View file

@ -5,6 +5,7 @@ module tally_filter_mesh
use hdf5
use constants
use dict_header, only: EMPTY
use error
use mesh_header, only: RegularMesh, meshes, n_meshes, mesh_dict
use hdf5_interface
@ -41,6 +42,7 @@ contains
integer :: i_mesh
integer :: id
integer :: n
integer :: val
n = node_word_count(node, "bins")
@ -51,8 +53,9 @@ contains
call get_node_value(node, "bins", id)
! Get pointer to mesh
if (mesh_dict % has_key(id)) then
i_mesh = mesh_dict % get_key(id)
val = mesh_dict % get(id)
if (val /= EMPTY) then
i_mesh = val
else
call fatal_error("Could not find mesh " // trim(to_str(id)) &
// " specified on filter.")

7
src/tallies/tally_filter_surface.F90
View file

@ -5,6 +5,7 @@ module tally_filter_surface
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error
use hdf5_interface
use surface_header
@ -84,12 +85,14 @@ contains
class(SurfaceFilter), intent(inout) :: this
integer :: i, id
integer :: val
! Convert ids to indices.
do i = 1, this % n_bins
id = this % surfaces(i)
if (surface_dict % has_key(id)) then
this % surfaces(i) = surface_dict % get_key(id)
val = surface_dict % get(id)
if (val /= EMPTY) then
this % surfaces(i) = val
else
call fatal_error("Could not find surface " // trim(to_str(id)) &
&// " specified on tally filter.")

16
src/tallies/tally_filter_universe.F90
View file

@ -5,6 +5,7 @@ module tally_filter_universe
use hdf5
use constants, only: ONE, MAX_LINE_LEN
use dict_header, only: EMPTY
use error, only: fatal_error
use hdf5_interface
use geometry_header
@ -54,12 +55,13 @@ contains
type(TallyFilterMatch), intent(inout) :: match
integer :: i
integer :: val
! Iterate over coordinate levels to see which universes match
do i = 1, p % n_coord
if (this % map % has_key(p % coord(i) % universe)) then
call match % bins % push_back(this % map % get_key(p % coord(i) &
% universe))
val = this % map % get(p % coord(i) % universe)
if (val /= EMPTY) then
call match % bins % push_back(val)
call match % weights % push_back(ONE)
end if
end do
@ -87,12 +89,14 @@ contains
class(UniverseFilter), intent(inout) :: this
integer :: i, id
integer :: val
! Convert ids to indices.
do i = 1, this % n_bins
id = this % universes(i)
if (universe_dict % has_key(id)) then
this % universes(i) = universe_dict % get_key(id)
val = universe_dict % get(id)
if (val /= EMPTY) then
this % universes(i) = val
else
call fatal_error("Could not find universe " // trim(to_str(id)) &
&// " specified on a tally filter.")
@ -101,7 +105,7 @@ contains
! Generate mapping from universe indices to filter bins.
do i = 1, this % n_bins
call this % map % add_key(this % universes(i), i)
call this % map % set(this % universes(i), i)
end do
end subroutine initialize_universe

10
src/tallies/tally_header.F90
View file

@ -455,8 +455,8 @@ contains
integer(C_INT) :: err
if (allocated(tallies)) then
if (tally_dict % has_key(id)) then
index = tally_dict % get_key(id)
if (tally_dict % has(id)) then
index = tally_dict % get(id)
err = 0
else
err = E_INVALID_ID
@ -592,7 +592,7 @@ contains
if (index >= 1 .and. index <= n_tallies) then
if (allocated(tallies(index) % obj)) then
tallies(index) % obj % id = id
call tally_dict % add_key(id, index)
call tally_dict % set(id, index)
err = 0
else
@ -633,8 +633,8 @@ contains
case ('total')
t % nuclide_bins(i) = -1
case default
if (nuclide_dict % has_key(nuclide_)) then
t % nuclide_bins(i) = nuclide_dict % get_key(nuclide_)
if (nuclide_dict % has(nuclide_)) then
t % nuclide_bins(i) = nuclide_dict % get(nuclide_)
else
err = E_DATA
call set_errmsg("Nuclide '" // trim(to_f_string(string)) // &

8
tests/test_enrichment/test_enrichment.py
View file

@ -21,7 +21,7 @@ if __name__ == '__main__':
sum_densities = 0.
for nuc in densities.keys():
assert nuc in ('U234', 'U235', 'U238')
assert nuc in ('U234', 'U235', 'U236', 'U238')
sum_densities += densities[nuc][1]
# Compute the weight percent U235
@ -30,4 +30,8 @@ if __name__ == '__main__':
# Compute the ratio of U234/U235
u234_to_u235 = densities['U234'][1] / densities['U235'][1]
assert np.isclose(u234_to_u235, 0.008, rtol=1.e-8)
assert np.isclose(u234_to_u235, 0.0089, rtol=1.e-8)
# Compute the ratio of U236/U235
u236_to_u235 = densities['U236'][1] / densities['U235'][1]
assert np.isclose(u236_to_u235, 0.0046, rtol=1.e-8)