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

Merge remote-tracking branch 'upstream/develop' into pyapi_lattice_outer

Browse source 
Conflicts:
	openmc/checkvalue.py
	openmc/filter.py
This commit is contained in:
Sterling Harper 2015-08-05 19:24:42 -06:00
commit 76e6133512
56 changed files with 1038 additions and 780 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/devguide/styleguide.rst
View file

@ -153,6 +153,8 @@ Avoid extraneous whitespace in the following situations:
Yes: if (variable == 2) then
No: if ( variable==2 ) then
Do not leave trailing whitespace at the end of a line.
------
Python
------

38
openmc/checkvalue.py
View file

@ -31,15 +31,15 @@ def check_type(name, value, expected_type, expected_iter_type=None):
"""
if not _isinstance(value, expected_type):
msg = 'Unable to set {0} to {1} which is not of type {2}'.format(
msg = 'Unable to set "{0}" to "{1}" which is not of type "{2}"'.format(
name, value, expected_type.__name__)
raise ValueError(msg)
if expected_iter_type:
for item in value:
if not _isinstance(item, expected_iter_type):
msg = 'Unable to set {0} to {1} since each item must be ' \
'of type {2}'.format(name, value,
msg = 'Unable to set "{0}" to "{1}" since each item must be ' \
'of type "{2}"'.format(name, value,
expected_iter_type.__name__)
raise ValueError(msg)
@ -89,7 +89,7 @@ def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1):
if _isinstance(current_item, expected_type):
# Is this deep enough?
if len(tree) < min_depth:
msg = 'Error setting {0}: The item at {1} does not meet the ' \
msg = 'Error setting "{0}": The item at {1} does not meet the '\
'minimum depth of {2}'.format(name, ind_str, min_depth)
raise ValueError(msg)
@ -139,16 +139,16 @@ def check_length(name, value, length_min, length_max=None):
if length_max is None:
if len(value) != length_min:
msg = 'Unable to set {0} to {1} since it must be of ' \
'length {2}'.format(name, value, length_min)
msg = 'Unable to set "{0}" to "{1}" since it must be of ' \
'length "{2}"'.format(name, value, length_min)
raise ValueError(msg)
elif not length_min <= len(value) <= length_max:
if length_min == length_max:
msg = 'Unable to set {0} to {1} since it must be of ' \
'length {2}'.format(name, value, length_min)
msg = 'Unable to set "{0}" to "{1}" since it must be of ' \
'length "{2}"'.format(name, value, length_min)
else:
msg = 'Unable to set {0} to {1} since it must have length ' \
'between {2} and {3}'.format(name, value, length_min,
msg = 'Unable to set "{0}" to "{1}" since it must have length ' \
'between "{2}" and "{3}"'.format(name, value, length_min,
length_max)
raise ValueError(msg)
@ -168,7 +168,7 @@ def check_value(name, value, accepted_values):
"""
if value not in accepted_values:
msg = 'Unable to set {0} to {1} since it is not in {2}'.format(
msg = 'Unable to set "{0}" to "{1}" since it is not in "{2}"'.format(
name, value, accepted_values)
raise ValueError(msg)
@ -190,13 +190,13 @@ def check_less_than(name, value, maximum, equality=False):
if equality:
if value > maximum:
msg = 'Unable to set {0} to {1} since it is greater than ' \
'{2}'.format(name, value, maximum)
msg = 'Unable to set "{0}" to "{1}" since it is greater than ' \
'"{2}"'.format(name, value, maximum)
raise ValueError(msg)
else:
if value >= maximum:
msg = 'Unable to set {0} to {1} since it is greater than ' \
'or equal to {2}'.format(name, value, maximum)
msg = 'Unable to set "{0}" to "{1}" since it is greater than ' \
'or equal to "{2}"'.format(name, value, maximum)
raise ValueError(msg)
def check_greater_than(name, value, minimum, equality=False):
@ -217,11 +217,11 @@ def check_greater_than(name, value, minimum, equality=False):
if equality:
if value < minimum:
msg = 'Unable to set {0} to {1} since it is less than ' \
'{2}'.format(name, value, minimum)
msg = 'Unable to set "{0}" to "{1}" since it is less than ' \
'"{2}"'.format(name, value, minimum)
raise ValueError(msg)
else:
if value <= minimum:
msg = 'Unable to set {0} to {1} since it is less than ' \
'or equal to {2}'.format(name, value, minimum)
msg = 'Unable to set "{0}" to "{1}" since it is less than ' \
'or equal to "{2}"'.format(name, value, minimum)
raise ValueError(msg)

2
openmc/executor.py
View file

@ -50,7 +50,7 @@ class Executor(object):
check_type("Executor's working directory", working_directory,
basestring)
if not os.path.isdir(working_directory):
msg = 'Unable to set Executor\'s working directory to {0} ' \
msg = 'Unable to set Executor\'s working directory to "{0}" ' \
'which does not exist'.format(working_directory)
raise ValueError(msg)

10
openmc/filter.py
View file

@ -142,19 +142,19 @@ class Filter(object):
elif self._type in ['energy', 'energyout']:
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \
msg = 'Unable to add bin edge "{0}" to a "{1}" Filter ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, self._type)
raise ValueError(msg)
elif edge < 0.:
msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \
msg = 'Unable to add bin edge "{0}" to a "{1}" Filter ' \
'since it is a negative value'.format(edge, self._type)
raise ValueError(msg)
# Check that bin edges are monotonically increasing
for index in range(len(bins)):
if index > 0 and bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a {1} Filter ' \
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, self._type)
raise ValueError(msg)
@ -201,7 +201,7 @@ class Filter(object):
def stride(self, stride):
check_type('filter stride', stride, Integral)
if stride < 0:
msg = 'Unable to set stride "{0}" for a {1} Filter since it is a ' \
msg = 'Unable to set stride "{0}" for a "{1}" Filter since it is a ' \
'negative value'.format(stride, self._type)
raise ValueError(msg)
@ -260,7 +260,7 @@ class Filter(object):
"""
if not self.can_merge(filter):
msg = 'Unable to merge {0} with {1} filters'.format(self._type, filter._type)
msg = 'Unable to merge "{0}" with "{1}" filters'.format(self._type, filter._type)
raise ValueError(msg)
# Create deep copy of filter to return as merged filter

4
openmc/geometry.py
View file

@ -34,8 +34,8 @@ class Geometry(object):
def root_universe(self, root_universe):
check_type('root universe', root_universe, openmc.Universe)
if root_universe._id != 0:
msg = 'Unable to add root Universe {0} to Geometry since ' \
'it has ID={1} instead of ' \
msg = 'Unable to add root Universe "{0}" to Geometry since ' \
'it has ID="{1}" instead of ' \
'ID=0'.format(root_universe, root_universe._id)
raise ValueError(msg)

56
openmc/material.py
View file

@ -26,7 +26,7 @@ def reset_auto_material_id():
# Units for density supported by OpenMC
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'at/b-cm', 'at/cm3', 'sum']
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum']
# Constant for density when not needed
NO_DENSITY = 99999.
@ -122,7 +122,7 @@ class Material(object):
else:
check_type('material ID', material_id, Integral)
if material_id in MATERIAL_IDS:
msg = 'Unable to set Material ID to {0} since a Material with ' \
msg = 'Unable to set Material ID to "{0}" since a Material with ' \
'this ID was already initialized'.format(material_id)
raise ValueError(msg)
check_greater_than('material ID', material_id, 0)
@ -132,7 +132,7 @@ class Material(object):
@name.setter
def name(self, name):
check_type('name for Material ID={0}'.format(self._id),
check_type('name for Material ID="{0}"'.format(self._id),
name, basestring)
self._name = name
@ -149,12 +149,12 @@ class Material(object):
"""
check_type('the density for Material ID={0}'.format(self._id),
check_type('the density for Material ID="{0}"'.format(self._id),
density, Real)
check_value('density units', units, DENSITY_UNITS)
if density == NO_DENSITY and units is not 'sum':
msg = 'Unable to set the density Material ID={0} ' \
msg = 'Unable to set the density Material ID="{0}" ' \
'because a density must be set when not using ' \
'sum unit'.format(self._id)
raise ValueError(msg)
@ -169,8 +169,8 @@ class Material(object):
'version of openmc')
if not isinstance(filename, basestring) and filename is not None:
msg = 'Unable to add OTF material file to Material ID={0} with a ' \
'non-string name {1}'.format(self._id, filename)
msg = 'Unable to add OTF material file to Material ID="{0}" with a ' \
'non-string name "{1}"'.format(self._id, filename)
raise ValueError(msg)
self._distrib_otf_file = filename
@ -198,18 +198,18 @@ class Material(object):
"""
if not isinstance(nuclide, (openmc.Nuclide, str)):
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'non-Nuclide value {1}'.format(self._id, nuclide)
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'non-Nuclide value "{1}"'.format(self._id, nuclide)
raise ValueError(msg)
elif not isinstance(percent, Real):
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'non-floating point value "{1}"'.format(self._id, percent)
raise ValueError(msg)
elif percent_type not in ['ao', 'wo', 'at/g-cm']:
msg = 'Unable to add a Nuclide to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
msg = 'Unable to add a Nuclide to Material ID="{0}" with a ' \
'percent type "{1}"'.format(self._id, percent_type)
raise ValueError(msg)
if isinstance(nuclide, openmc.Nuclide):
@ -232,7 +232,7 @@ class Material(object):
"""
if not isinstance(nuclide, openmc.Nuclide):
msg = 'Unable to remove a Nuclide {0} in Material ID={1} ' \
msg = 'Unable to remove a Nuclide "{0}" in Material ID="{1}" ' \
'since it is not a Nuclide'.format(self._id, nuclide)
raise ValueError(msg)
@ -255,18 +255,18 @@ class Material(object):
"""
if not isinstance(element, openmc.Element):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-Element value {1}'.format(self._id, element)
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-Element value "{1}"'.format(self._id, element)
raise ValueError(msg)
if not isinstance(percent, Real):
msg = 'Unable to add an Element to Material ID={0} with a ' \
'non-floating point value {1}'.format(self._id, percent)
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'non-floating point value "{1}"'.format(self._id, percent)
raise ValueError(msg)
if percent_type not in ['ao', 'wo']:
msg = 'Unable to add an Element to Material ID={0} with a ' \
'percent type {1}'.format(self._id, percent_type)
msg = 'Unable to add an Element to Material ID="{0}" with a ' \
'percent type "{1}"'.format(self._id, percent_type)
raise ValueError(msg)
# Copy this Element to separate it from same Element in other Materials
@ -301,13 +301,13 @@ class Material(object):
"""
if not isinstance(name, basestring):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string table name {1}'.format(self._id, name)
msg = 'Unable to add an S(a,b) table to Material ID="{0}" with a ' \
'non-string table name "{1}"'.format(self._id, name)
raise ValueError(msg)
if not isinstance(xs, basestring):
msg = 'Unable to add an S(a,b) table to Material ID={0} with a ' \
'non-string cross-section identifier {1}'.format(self._id, xs)
msg = 'Unable to add an S(a,b) table to Material ID="{0}" with a ' \
'non-string cross-section identifier "{1}"'.format(self._id, xs)
raise ValueError(msg)
self._sab.append((name, xs))
@ -332,7 +332,7 @@ class Material(object):
return nuclides
def _repr__(self):
def __repr__(self):
string = 'Material\n'
string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id)
string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name)
@ -522,7 +522,7 @@ class MaterialsFile(object):
"""
if not isinstance(material, Material):
msg = 'Unable to add a non-Material {0} to the ' \
msg = 'Unable to add a non-Material "{0}" to the ' \
'MaterialsFile'.format(material)
raise ValueError(msg)
@ -539,7 +539,7 @@ class MaterialsFile(object):
"""
if not isinstance(materials, Iterable):
msg = 'Unable to create OpenMC materials.xml file from {0} which ' \
msg = 'Unable to create OpenMC materials.xml file from "{0}" which ' \
'is not iterable'.format(materials)
raise ValueError(msg)
@ -557,7 +557,7 @@ class MaterialsFile(object):
"""
if not isinstance(material, Material):
msg = 'Unable to remove a non-Material {0} from the ' \
msg = 'Unable to remove a non-Material "{0}" from the ' \
'MaterialsFile'.format(material)
raise ValueError(msg)

6
openmc/mesh.py
View file

@ -148,14 +148,14 @@ class Mesh(object):
@name.setter
def name(self, name):
check_type('name for mesh ID={0}'.format(self._id), name, basestring)
check_type('name for mesh ID="{0}"'.format(self._id), name, basestring)
self._name = name
@type.setter
def type(self, meshtype):
check_type('type for mesh ID={0}'.format(self._id),
check_type('type for mesh ID="{0}"'.format(self._id),
meshtype, basestring)
check_value('type for mesh ID={0}'.format(self._id),
check_value('type for mesh ID="{0}"'.format(self._id),
meshtype, ['rectangular', 'hexagonal'])
self._type = meshtype

40
openmc/opencg_compatible.py
View file

@ -78,7 +78,7 @@ def get_opencg_material(openmc_material):
"""
if not isinstance(openmc_material, openmc.Material):
msg = 'Unable to create an OpenCG Material from {0} ' \
msg = 'Unable to create an OpenCG Material from "{0}" ' \
'which is not an OpenMC Material'.format(openmc_material)
raise ValueError(msg)
@ -118,7 +118,7 @@ def get_openmc_material(opencg_material):
"""
if not isinstance(opencg_material, opencg.Material):
msg = 'Unable to create an OpenMC Material from {0} ' \
msg = 'Unable to create an OpenMC Material from "{0}" ' \
'which is not an OpenCG Material'.format(opencg_material)
raise ValueError(msg)
@ -165,7 +165,7 @@ def is_opencg_surface_compatible(opencg_surface):
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to check if OpenCG Surface is compatible' \
'since {0} is not a Surface'.format(opencg_surface)
'since "{0}" is not a Surface'.format(opencg_surface)
raise ValueError(msg)
if opencg_surface._type in ['x-squareprism',
@ -191,7 +191,7 @@ def get_opencg_surface(openmc_surface):
"""
if not isinstance(openmc_surface, openmc.Surface):
msg = 'Unable to create an OpenCG Surface from {0} ' \
msg = 'Unable to create an OpenCG Surface from "{0}" ' \
'which is not an OpenMC Surface'.format(openmc_surface)
raise ValueError(msg)
@ -277,7 +277,7 @@ def get_openmc_surface(opencg_surface):
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
msg = 'Unable to create an OpenMC Surface from "{0}" which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
@ -335,7 +335,7 @@ def get_openmc_surface(opencg_surface):
else:
msg = 'Unable to create an OpenMC Surface from an OpenCG ' \
'Surface of type {0} since it is not a compatible ' \
'Surface of type "{0}" since it is not a compatible ' \
'Surface type in OpenMC'.format(opencg_surface._type)
raise ValueError(msg)
@ -368,7 +368,7 @@ def get_compatible_opencg_surfaces(opencg_surface):
"""
if not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create an OpenMC Surface from {0} which ' \
msg = 'Unable to create an OpenMC Surface from "{0}" which ' \
'is not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
@ -421,7 +421,7 @@ def get_compatible_opencg_surfaces(opencg_surface):
else:
msg = 'Unable to create a compatible OpenMC Surface an OpenCG ' \
'Surface of type {0} since it already a compatible ' \
'Surface of type "{0}" since it already a compatible ' \
'Surface type in OpenMC'.format(opencg_surface._type)
raise ValueError(msg)
@ -450,7 +450,7 @@ def get_opencg_cell(openmc_cell):
"""
if not isinstance(openmc_cell, openmc.Cell):
msg = 'Unable to create an OpenCG Cell from {0} which ' \
msg = 'Unable to create an OpenCG Cell from "{0}" which ' \
'is not an OpenMC Cell'.format(openmc_cell)
raise ValueError(msg)
@ -518,17 +518,17 @@ def get_compatible_opencg_cells(opencg_cell, opencg_surface, halfspace):
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create compatible OpenMC Cell from {0} which ' \
msg = 'Unable to create compatible OpenMC Cell from "{0}" which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
elif not isinstance(opencg_surface, opencg.Surface):
msg = 'Unable to create compatible OpenMC Cell since {0} is ' \
msg = 'Unable to create compatible OpenMC Cell since "{0}" is ' \
'not an OpenCG Surface'.format(opencg_surface)
raise ValueError(msg)
elif halfspace not in [-1, +1]:
msg = 'Unable to create compatible Cell since {0}' \
msg = 'Unable to create compatible Cell since "{0}"' \
'is not a +/-1 halfspace'.format(halfspace)
raise ValueError(msg)
@ -626,7 +626,7 @@ def make_opencg_cells_compatible(opencg_universe):
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to make compatible OpenCG Cells for {0} which ' \
msg = 'Unable to make compatible OpenCG Cells for "{0}" which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
@ -685,7 +685,7 @@ def get_openmc_cell(opencg_cell):
"""
if not isinstance(opencg_cell, opencg.Cell):
msg = 'Unable to create an OpenMC Cell from {0} which ' \
msg = 'Unable to create an OpenMC Cell from "{0}" which ' \
'is not an OpenCG Cell'.format(opencg_cell)
raise ValueError(msg)
@ -749,7 +749,7 @@ def get_opencg_universe(openmc_universe):
"""
if not isinstance(openmc_universe, openmc.Universe):
msg = 'Unable to create an OpenCG Universe from {0} which ' \
msg = 'Unable to create an OpenCG Universe from "{0}" which ' \
'is not an OpenMC Universe'.format(openmc_universe)
raise ValueError(msg)
@ -796,7 +796,7 @@ def get_openmc_universe(opencg_universe):
"""
if not isinstance(opencg_universe, opencg.Universe):
msg = 'Unable to create an OpenMC Universe from {0} which ' \
msg = 'Unable to create an OpenMC Universe from "{0}" which ' \
'is not an OpenCG Universe'.format(opencg_universe)
raise ValueError(msg)
@ -846,7 +846,7 @@ def get_opencg_lattice(openmc_lattice):
"""
if not isinstance(openmc_lattice, openmc.Lattice):
msg = 'Unable to create an OpenCG Lattice from {0} which ' \
msg = 'Unable to create an OpenCG Lattice from "{0}" which ' \
'is not an OpenMC Lattice'.format(openmc_lattice)
raise ValueError(msg)
@ -926,7 +926,7 @@ def get_openmc_lattice(opencg_lattice):
"""
if not isinstance(opencg_lattice, opencg.Lattice):
msg = 'Unable to create an OpenMC Lattice from {0} which ' \
msg = 'Unable to create an OpenMC Lattice from "{0}" which ' \
'is not an OpenCG Lattice'.format(opencg_lattice)
raise ValueError(msg)
@ -997,7 +997,7 @@ def get_opencg_geometry(openmc_geometry):
"""
if not isinstance(openmc_geometry, openmc.Geometry):
msg = 'Unable to get OpenCG geometry from {0} which is ' \
msg = 'Unable to get OpenCG geometry from "{0}" which is ' \
'not an OpenMC Geometry object'.format(openmc_geometry)
raise ValueError(msg)
@ -1037,7 +1037,7 @@ def get_openmc_geometry(opencg_geometry):
"""
if not isinstance(opencg_geometry, opencg.Geometry):
msg = 'Unable to get OpenMC geometry from {0} which is ' \
msg = 'Unable to get OpenMC geometry from "{0}" which is ' \
'not an OpenCG Geometry object'.format(opencg_geometry)
raise ValueError(msg)

12
openmc/plots.py
View file

@ -210,18 +210,18 @@ class Plot(object):
for key in col_spec:
if key < 0:
msg = 'Unable to create Plot ID={0} with col_spec ID {1} ' \
msg = 'Unable to create Plot ID="{0}" with col_spec ID "{1}" ' \
'which is less than 0'.format(self._id, key)
raise ValueError(msg)
elif not isinstance(col_spec[key], Iterable):
msg = 'Unable to create Plot ID={0} with col_spec RGB values' \
' {1} which is not iterable'.format(self._id, col_spec[key])
msg = 'Unable to create Plot ID="{0}" with col_spec RGB values' \
' "{1}" which is not iterable'.format(self._id, col_spec[key])
raise ValueError(msg)
elif len(col_spec[key]) != 3:
msg = 'Unable to create Plot ID={0} with col_spec RGB ' \
'values of length {1} since 3 values must be ' \
msg = 'Unable to create Plot ID="{0}" with col_spec RGB ' \
'values of length "{1}" since 3 values must be ' \
'input'.format(self._id, len(col_spec[key]))
raise ValueError(msg)
@ -332,7 +332,7 @@ class PlotsFile(object):
"""
if not isinstance(plot, Plot):
msg = 'Unable to add a non-Plot {0} to the PlotsFile'.format(plot)
msg = 'Unable to add a non-Plot "{0}" to the PlotsFile'.format(plot)
raise ValueError(msg)
self._plots.append(plot)

22
openmc/settings.py
View file

@ -426,28 +426,28 @@ class SettingsFile(object):
@keff_trigger.setter
def keff_trigger(self, keff_trigger):
if not isinstance(keff_trigger, dict):
msg = 'Unable to set a trigger on keff from {0} which ' \
msg = 'Unable to set a trigger on keff from "{0}" which ' \
'is not a Python dictionary'.format(keff_trigger)
raise ValueError(msg)
elif 'type' not in keff_trigger:
msg = 'Unable to set a trigger on keff from {0} which ' \
msg = 'Unable to set a trigger on keff from "{0}" which ' \
'does not have a "type" key'.format(keff_trigger)
raise ValueError(msg)
elif keff_trigger['type'] not in ['variance', 'std_dev', 'rel_err']:
msg = 'Unable to set a trigger on keff with ' \
'type {0}'.format(keff_trigger['type'])
'type "{0}"'.format(keff_trigger['type'])
raise ValueError(msg)
elif 'threshold' not in keff_trigger:
msg = 'Unable to set a trigger on keff from {0} which ' \
msg = 'Unable to set a trigger on keff from "{0}" which ' \
'does not have a "threshold" key'.format(keff_trigger)
raise ValueError(msg)
elif not isinstance(keff_trigger['threshold'], Real):
msg = 'Unable to set a trigger on keff with ' \
'threshold {0}'.format(keff_trigger['threshold'])
'threshold "{0}"'.format(keff_trigger['threshold'])
raise ValueError(msg)
self._keff_trigger = keff_trigger
@ -574,20 +574,20 @@ class SettingsFile(object):
@output.setter
def output(self, output):
if not isinstance(output, dict):
msg = 'Unable to set output to {0} which is not a Python ' \
msg = 'Unable to set output to "{0}" which is not a Python ' \
'dictionary of string keys and boolean values'.format(output)
raise ValueError(msg)
for element in output:
keys = ['summary', 'cross_sections', 'tallies', 'distribmats']
if element not in keys:
msg = 'Unable to set output to {0} which is unsupported by ' \
msg = 'Unable to set output to "{0}" which is unsupported by ' \
'OpenMC'.format(element)
raise ValueError(msg)
if not isinstance(output[element], (bool, np.bool)):
msg = 'Unable to set output for {0} to a non-boolean ' \
'value {1}'.format(element, output[element])
msg = 'Unable to set output for "{0}" to a non-boolean ' \
'value "{1}"'.format(element, output[element])
raise ValueError(msg)
self._output = output
@ -753,7 +753,7 @@ class SettingsFile(object):
def track(self, track):
check_type('track', track, Iterable, Integral)
if len(track) % 3 != 0:
msg = 'Unable to set the track to {0} since its length is ' \
msg = 'Unable to set the track to "{0}" since its length is ' \
'not a multiple of 3'.format(track)
raise ValueError(msg)
for t in zip(track[::3], track[1::3], track[2::3]):
@ -832,7 +832,7 @@ class SettingsFile(object):
len_nodemap = np.prod(self._dd_mesh_dimension)
if len(nodemap) < len_nodemap or len(nodemap) > len_nodemap:
msg = 'Unable to set DD nodemap with length {0} which ' \
msg = 'Unable to set DD nodemap with length "{0}" which ' \
'does not have the same dimensionality as the domain ' \
'mesh'.format(len(nodemap))
raise ValueError(msg)

4
openmc/statepoint.py
View file

@ -372,7 +372,7 @@ class StatePoint(object):
path='{0}{1}/n_bins'.format(subbase, j))[0]
if n_bins <= 0:
msg = 'Unable to create Filter {0} for Tally ID={2} ' \
msg = 'Unable to create Filter "{0}" for Tally ID="{1}" ' \
'since no bins were specified'.format(j, tally_key)
raise ValueError(msg)
@ -748,7 +748,7 @@ class StatePoint(object):
"""
if not isinstance(summary, openmc.summary.Summary):
msg = 'Unable to link statepoint with {0} which ' \
msg = 'Unable to link statepoint with "{0}" which ' \
'is not a Summary object'.format(summary)
raise ValueError(msg)

40
openmc/tallies.py
View file

@ -297,7 +297,7 @@ class Tally(object):
"""
if not isinstance(trigger, Trigger):
msg = 'Unable to add a tally trigger for Tally ID={0} to ' \
msg = 'Unable to add a tally trigger for Tally ID="{0}" to ' \
'since "{1}" is not a Trigger'.format(self.id, trigger)
raise ValueError(msg)
@ -330,7 +330,7 @@ class Tally(object):
"""
if not isinstance(filter, Filter):
msg = 'Unable to add Filter "{0}" to Tally ID={1} since it is ' \
msg = 'Unable to add Filter "{0}" to Tally ID="{1}" since it is ' \
'not a Filter object'.format(filter, self.id)
raise ValueError(msg)
@ -359,7 +359,7 @@ class Tally(object):
"""
if not isinstance(score, basestring):
msg = 'Unable to add score "{0}" to Tally ID={1} since it is ' \
msg = 'Unable to add score "{0}" to Tally ID="{1}" since it is ' \
'not a string'.format(score, self.id)
raise ValueError(msg)
@ -405,7 +405,7 @@ class Tally(object):
"""
if score not in self.scores:
msg = 'Unable to remove score "{0}" from Tally ID={1} since the ' \
msg = 'Unable to remove score "{0}" from Tally ID="{1}" since the ' \
'Tally does not contain this score'.format(score, self.id)
ValueError(msg)
@ -422,7 +422,7 @@ class Tally(object):
"""
if filter not in self.filters:
msg = 'Unable to remove filter "{0}" from Tally ID={1} since the ' \
msg = 'Unable to remove filter "{0}" from Tally ID="{1}" since the ' \
'Tally does not contain this filter'.format(filter, self.id)
ValueError(msg)
@ -439,7 +439,7 @@ class Tally(object):
"""
if nuclide not in self.nuclides:
msg = 'Unable to remove nuclide "{0}" from Tally ID={1} since the ' \
msg = 'Unable to remove nuclide "{0}" from Tally ID="{1}" since the ' \
'Tally does not contain this nuclide'.format(nuclide, self.id)
ValueError(msg)
@ -555,7 +555,7 @@ class Tally(object):
"""
if not self.can_merge(tally):
msg = 'Unable to merge tally ID={0} with {1}'.format(tally.id, self.id)
msg = 'Unable to merge tally ID="{0}" with "{1}"'.format(tally.id, self.id)
raise ValueError(msg)
# Create deep copy of tally to return as merged tally
@ -627,7 +627,7 @@ class Tally(object):
# Scores
if len(self.scores) == 0:
msg = 'Unable to get XML for Tally ID={0} since it does not ' \
msg = 'Unable to get XML for Tally ID="{0}" since it does not ' \
'contain any scores'.format(self.id)
raise ValueError(msg)
@ -682,7 +682,7 @@ class Tally(object):
# If we did not find the Filter, throw an Exception
if filter is None:
msg = 'Unable to find filter type "{0}" in ' \
'Tally ID={1}'.format(filter_type, self.id)
'Tally ID="{1}"'.format(filter_type, self.id)
raise ValueError(msg)
return filter
@ -849,7 +849,7 @@ class Tally(object):
# Ensure that StatePoint.read_results() was called first
if self._sum is None or self._sum_sq is None:
msg = 'The Tally ID={0} has no data to return. Call the ' \
msg = 'The Tally ID="{0}" has no data to return. Call the ' \
'StatePoint.read_results() routine before using ' \
'Tally.get_values(...)'.format(self.id)
raise ValueError(msg)
@ -944,7 +944,7 @@ class Tally(object):
elif value == 'sum_sq':
data = self.sum_sq[indices]
else:
msg = 'Unable to return results from Tally ID={0} since the ' \
msg = 'Unable to return results from Tally ID="{0}" since the ' \
'the requested value "{1}" is not \'mean\', \'std_dev\', ' \
'\rel_err\', \'sum\', or \'sum_sq\''.format(self.id, value)
raise LookupError(msg)
@ -996,14 +996,14 @@ class Tally(object):
# Ensure that StatePoint.read_results() was called first
if self._sum is None or self._sum_sq is None:
msg = 'The Tally ID={0} has no data to return. Call the ' \
msg = 'The Tally ID="{0}" has no data to return. Call the ' \
'StatePoint.read_results() routine before using ' \
'Tally.get_pandas_dataframe(...)'.format(self.id)
raise KeyError(msg)
# If using Summary, ensure StatePoint.link_with_summary(...) was called
if summary and not self.with_summary:
msg = 'The Tally ID={0} has not been linked with the Summary. ' \
msg = 'The Tally ID="{0}" has not been linked with the Summary. ' \
'Call the StatePoint.link_with_summary(...) routine ' \
'before using Tally.get_pandas_dataframe(...) with ' \
'Summary info'.format(self.id)
@ -1292,30 +1292,30 @@ class Tally(object):
# Ensure that StatePoint.read_results() was called first
if self._sum is None or self._sum_sq is None:
msg = 'The Tally ID={0} has no data to export. Call the ' \
msg = 'The Tally ID="{0}" has no data to export. Call the ' \
'StatePoint.read_results() routine before using ' \
'Tally.export_results(...)'.format(self.id)
raise KeyError(msg)
if not isinstance(filename, basestring):
msg = 'Unable to export the results for Tally ID={0} to ' \
msg = 'Unable to export the results for Tally ID="{0}" to ' \
'filename="{1}" since it is not a ' \
'string'.format(self.id, filename)
raise ValueError(msg)
elif not isinstance(directory, basestring):
msg = 'Unable to export the results for Tally ID={0} to ' \
msg = 'Unable to export the results for Tally ID="{0}" to ' \
'directory="{1}" since it is not a ' \
'string'.format(self.id, directory)
raise ValueError(msg)
elif format not in ['hdf5', 'pkl', 'csv']:
msg = 'Unable to export the results for Tally ID={0} to format ' \
msg = 'Unable to export the results for Tally ID="{0}" to format ' \
'"{1}" since it is not supported'.format(self.id, format)
raise ValueError(msg)
elif not isinstance(append, bool):
msg = 'Unable to export the results for Tally ID={0} since the ' \
msg = 'Unable to export the results for Tally ID="{0}" since the ' \
'append parameter is not True/False'.format(self.id, append)
raise ValueError(msg)
@ -1437,7 +1437,7 @@ class TalliesFile(object):
"""
if not isinstance(tally, Tally):
msg = 'Unable to add a non-Tally {0} to the TalliesFile'.format(tally)
msg = 'Unable to add a non-Tally "{0}" to the TalliesFile'.format(tally)
raise ValueError(msg)
if merge:
@ -1508,7 +1508,7 @@ class TalliesFile(object):
"""
if not isinstance(mesh, Mesh):
msg = 'Unable to add a non-Mesh {0} to the TalliesFile'.format(mesh)
msg = 'Unable to add a non-Mesh "{0}" to the TalliesFile'.format(mesh)
raise ValueError(msg)
self._meshes.append(mesh)

22
openmc/universe.py
View file

@ -127,8 +127,8 @@ class Cell(object):
if fill.strip().lower() == 'void':
self._type = 'void'
else:
msg = 'Unable to set Cell ID={0} to use a non-Material or ' \
'Universe fill {1}'.format(self._id, fill)
msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
'Universe fill "{1}"'.format(self._id, fill)
raise ValueError(msg)
elif isinstance(fill, openmc.Material):
@ -141,8 +141,8 @@ class Cell(object):
self._type = 'lattice'
else:
msg = 'Unable to set Cell ID={0} to use a non-Material or ' \
'Universe fill {1}'.format(self._id, fill)
msg = 'Unable to set Cell ID="{0}" to use a non-Material or ' \
'Universe fill "{1}"'.format(self._id, fill)
raise ValueError(msg)
self._fill = fill
@ -178,13 +178,13 @@ class Cell(object):
"""
if not isinstance(surface, openmc.Surface):
msg = 'Unable to add Surface {0} to Cell ID={1} since it is ' \
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
raise ValueError(msg)
if halfspace not in [-1, +1]:
msg = 'Unable to add Surface {0} to Cell ID={1} with halfspace ' \
'{2} since it is not +/-1'.format(surface, self._id, halfspace)
msg = 'Unable to add Surface "{0}" to Cell ID="{1}" with halfspace ' \
'"{2}" since it is not +/-1'.format(surface, self._id, halfspace)
raise ValueError(msg)
# If the Cell does not already contain the Surface, add it
@ -202,7 +202,7 @@ class Cell(object):
"""
if not isinstance(surface, openmc.Surface):
msg = 'Unable to remove Surface {0} from Cell ID={1} since it is ' \
msg = 'Unable to remove Surface "{0}" from Cell ID="{1}" since it is ' \
'not a Surface object'.format(surface, self._id)
raise ValueError(msg)
@ -453,7 +453,7 @@ class Universe(object):
"""
if not isinstance(cell, Cell):
msg = 'Unable to add a Cell to Universe ID={0} since {1} is not ' \
msg = 'Unable to add a Cell to Universe ID="{0}" since "{1}" is not ' \
'a Cell'.format(self._id, cell)
raise ValueError(msg)
@ -473,7 +473,7 @@ class Universe(object):
"""
if not isinstance(cells, Iterable):
msg = 'Unable to add Cells to Universe ID={0} since {1} is not ' \
msg = 'Unable to add Cells to Universe ID="{0}" since "{1}" is not ' \
'iterable'.format(self._id, cells)
raise ValueError(msg)
@ -491,7 +491,7 @@ class Universe(object):
"""
if not isinstance(cell, Cell):
msg = 'Unable to remove a Cell from Universe ID={0} since {1} is ' \
msg = 'Unable to remove a Cell from Universe ID="{0}" since "{1}" is ' \
'not a Cell'.format(self._id, cell)
raise ValueError(msg)

27
src/ace.F90
View file

@ -88,16 +88,16 @@ contains
nuclides(i_nuclide) % resonant = .true.
nuclides(i_nuclide) % name_0K = nuclides_0K(n) % name_0K
nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % &
& name_0K)
& name_0K)
nuclides(i_nuclide) % scheme = nuclides_0K(n) % scheme
nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % &
& scheme)
& scheme)
nuclides(i_nuclide) % E_min = nuclides_0K(n) % E_min
nuclides(i_nuclide) % E_max = nuclides_0K(n) % E_max
if (.not. already_read % contains(nuclides(i_nuclide) % &
& name_0K)) then
& name_0K)) then
i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % &
& name_0K)
& name_0K)
call read_ace_table(i_nuclide, i_listing)
end if
exit
@ -446,9 +446,10 @@ contains
if (nuc % elastic_0K(i) < ZERO) nuc % elastic_0K(i) = ZERO
! build xs cdf
xs_cdf_sum = xs_cdf_sum + (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) &
& + sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO &
& * (nuc % energy_0K(i+1) - nuc % energy_0K(i))
xs_cdf_sum = xs_cdf_sum &
+ (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) &
+ sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO &
* (nuc % energy_0K(i+1) - nuc % energy_0K(i))
nuc % xs_cdf(i) = xs_cdf_sum
end do
@ -752,31 +753,31 @@ contains
! Set flag and allocate space for Tab1 to store yield
rxn % multiplicity_with_E = .true.
allocate(rxn % multiplicity_E)
XSS_index = JXS(11) + rxn % multiplicity - 101
NR = nint(XSS(XSS_index))
rxn % multiplicity_E % n_regions = NR
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(rxn % multiplicity_E % nbt(NR))
allocate(rxn % multiplicity_E % int(NR))
end if
! read ENDF interpolation parameters
XSS_index = XSS_index + 1
if (NR > 0) then
rxn % multiplicity_E % nbt = get_int(NR)
rxn % multiplicity_E % int = get_int(NR)
end if
! allocate space for yield data
XSS_index = XSS_index + 2*NR
NE = nint(XSS(XSS_index))
rxn % multiplicity_E % n_pairs = NE
allocate(rxn % multiplicity_E % x(NE))
allocate(rxn % multiplicity_E % y(NE))
! read yield data
XSS_index = XSS_index + 1
rxn % multiplicity_E % x = get_real(NE)
@ -1480,7 +1481,7 @@ contains
table % inelastic_data(i) % e_out_pdf(j) = XSS(XSS_index + 2)
table % inelastic_data(i) % e_out_cdf(j) = XSS(XSS_index + 3)
table % inelastic_data(i) % mu(:, j) = &
XSS(XSS_index + 4: XSS_index + 4 + NMU - 1)
XSS(XSS_index + 4: XSS_index + 4 + NMU - 1)
XSS_index = XSS_index + 4 + NMU - 1
end do
end do

30
src/cmfd_data.F90
View file

@ -104,23 +104,23 @@ contains
cmfd % keff_bal = ZERO
! Begin loop around tallies
TAL: do ital = 1, n_cmfd_tallies
! Begin loop around tallies
TAL: do ital = 1, n_cmfd_tallies
! Associate tallies and mesh
t => cmfd_tallies(ital)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
m => meshes(i_mesh)
! Associate tallies and mesh
t => cmfd_tallies(ital)
i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1)
m => meshes(i_mesh)
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
i_filter_eout = t % find_filter(FILTER_ENERGYOUT)
i_filter_surf = t % find_filter(FILTER_SURFACE)
i_filter_mesh = t % find_filter(FILTER_MESH)
i_filter_ein = t % find_filter(FILTER_ENERGYIN)
i_filter_eout = t % find_filter(FILTER_ENERGYOUT)
i_filter_surf = t % find_filter(FILTER_SURFACE)
! Begin loop around space
ZLOOP: do k = 1,nz
! Begin loop around space
ZLOOP: do k = 1,nz
YLOOP: do j = 1,ny
YLOOP: do j = 1,ny
XLOOP: do i = 1,nx
@ -600,7 +600,7 @@ contains
dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + &
cell_hxyz(xyz_idx)*neig_dc)
end if
end if
end if
@ -797,7 +797,7 @@ contains
! Calculate albedo
if ((shift_idx == 1 .and. current(2*l ) < 1.0e-10_8) .or. &
(shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then
(shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then
albedo = ONE
else
albedo = (current(2*l-1)/current(2*l))**(shift_idx)

2
src/cmfd_execute.F90
View file

@ -121,7 +121,7 @@ contains
! Allocate cmfd source if not already allocated and allocate buffer
if (.not. allocated(cmfd % cmfd_src)) &
allocate(cmfd % cmfd_src(ng,nx,ny,nz))
allocate(cmfd % cmfd_src(ng,nx,ny,nz))
! Reset cmfd source to 0
cmfd % cmfd_src = ZERO

10
src/cmfd_header.F90
View file

@ -1,4 +1,4 @@
module cmfd_header
module cmfd_header
use constants, only: CMFD_NOACCEL, ZERO, ONE
@ -17,7 +17,7 @@ module cmfd_header
! Core overlay map
integer, allocatable :: coremap(:,:,:)
integer, allocatable :: indexmap(:,:)
integer :: mat_dim = CMFD_NOACCEL
integer :: mat_dim = CMFD_NOACCEL
! Energy grid
real(8), allocatable :: egrid(:)
@ -51,7 +51,7 @@ module cmfd_header
! Source sites in each mesh box
real(8), allocatable :: sourcecounts(:,:,:,:)
! Weight adjustment factors
! Weight adjustment factors
real(8), allocatable :: weightfactors(:,:,:,:)
! Eigenvector/eigenvalue from cmfd run
@ -118,7 +118,7 @@ contains
if (.not. allocated(this % nfissxs)) allocate(this % nfissxs(ng,ng,nx,ny,nz))
if (.not. allocated(this % diffcof)) allocate(this % diffcof(ng,nx,ny,nz))
! Allocate dtilde and dhat
! Allocate dtilde and dhat
if (.not. allocated(this % dtilde)) allocate(this % dtilde(6,ng,nx,ny,nz))
if (.not. allocated(this % dhat)) allocate(this % dhat(6,ng,nx,ny,nz))
@ -167,7 +167,7 @@ contains
end subroutine allocate_cmfd
!===============================================================================
! DEALLOCATE_CMFD frees all memory of cmfd type
! DEALLOCATE_CMFD frees all memory of cmfd type
!===============================================================================
subroutine deallocate_cmfd(this)

18
src/cmfd_input.F90
View file

@ -120,7 +120,7 @@ contains
allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), &
cmfd % indices(3)))
if (get_arraysize_integer(node_mesh, "map") /= &
product(cmfd % indices(1:3))) then
product(cmfd % indices(1:3))) then
call fatal_error('CMFD coremap not to correct dimensions')
end if
allocate(iarray(get_arraysize_integer(node_mesh, "map")))
@ -156,7 +156,7 @@ contains
call get_node_value(doc, "dhat_reset", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
dhat_reset = .true.
dhat_reset = .true.
end if
! Set monitoring
@ -180,7 +180,7 @@ contains
call get_node_value(doc, "run_adjoint", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
cmfd_run_adjoint = .true.
cmfd_run_adjoint = .true.
end if
! Batch to begin cmfd
@ -289,7 +289,7 @@ contains
! Determine number of dimensions for mesh
n = get_arraysize_integer(node_mesh, "dimension")
if (n /= 2 .and. n /= 3) then
call fatal_error("Mesh must be two or three dimensions.")
call fatal_error("Mesh must be two or three dimensions.")
end if
m % n_dimension = n
@ -318,14 +318,14 @@ contains
! Make sure both upper-right or width were specified
if (check_for_node(node_mesh, "upper_right") .and. &
check_for_node(node_mesh, "width")) then
check_for_node(node_mesh, "width")) then
call fatal_error("Cannot specify both <upper_right> and <width> on a &
&tally mesh.")
end if
! Make sure either upper-right or width was specified
if (.not.check_for_node(node_mesh, "upper_right") .and. &
.not.check_for_node(node_mesh, "width")) then
.not.check_for_node(node_mesh, "width")) then
call fatal_error("Must specify either <upper_right> and <width> on a &
&tally mesh.")
end if
@ -333,7 +333,7 @@ contains
if (check_for_node(node_mesh, "width")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "width") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <width> must be the same as the &
&number of entries on <lower_left>.")
end if
@ -351,7 +351,7 @@ contains
elseif (check_for_node(node_mesh, "upper_right")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "upper_right") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <upper_right> must be the same &
&as the number of entries on <lower_left>.")
end if
@ -548,7 +548,7 @@ contains
! Deallocate filter bins
deallocate(filters(1) % int_bins)
if (check_for_node(node_mesh, "energy")) &
deallocate(filters(2) % real_bins)
deallocate(filters(2) % real_bins)
end do

20
src/cmfd_loss_operator.F90
View file

@ -54,7 +54,7 @@ contains
n_e = 4*(nx - 2) + 4*(ny - 2) + 4*(nz - 2) ! define # of edges
n_s = 2*(nx - 2)*(ny - 2) + 2*(nx - 2)*(nz - 2) &
+ 2*(ny - 2)*(nz - 2) ! define # of sides
n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors
n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors
nz_c = ng*n_c*(4 + ng - 1) ! define # nonzero corners
nz_e = ng*n_e*(5 + ng - 1) ! define # nonzero edges
nz_s = ng*n_s*(6 + ng - 1) ! define # nonzero sides
@ -131,7 +131,7 @@ contains
! Check for global boundary
if (bound(l) /= nxyz(xyz_idx,dir_idx)) then
! Check for coremap
! Check for coremap
if (cmfd_coremap) then
! Check for neighbor that is non-acceleartred
@ -139,11 +139,11 @@ contains
CMFD_NOACCEL) then
! Get neighbor matrix index
call indices_to_matrix(g,neig_idx(1), neig_idx(2), &
call indices_to_matrix(g,neig_idx(1), neig_idx(2), &
neig_idx(3), neig_mat_idx, ng, nx, ny)
! Record nonzero
nnz = nnz + 1
nnz = nnz + 1
end if
@ -218,11 +218,11 @@ contains
real(8) :: jn ! direction dependent leakage coeff to neig
real(8) :: jo(6) ! leakage coeff in front of cell flux
real(8) :: jnet ! net leakage from jo
real(8) :: val ! temporary variable before saving to
real(8) :: val ! temporary variable before saving to
! Check for adjoint
adjoint_calc = .false.
if (present(adjoint)) adjoint_calc = adjoint
if (present(adjoint)) adjoint_calc = adjoint
! Get maximum number of cells in each direction
nx = cmfd%indices(1)
@ -257,11 +257,11 @@ contains
dhat = ZERO
end if
! Create boundary vector
! Create boundary vector
bound = (/i,i,j,j,k,k/)
! Begin loop over leakages
! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z
! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z
LEAK: do l = 1,6
! Define (x,y,z) and (-,+) indices
@ -350,7 +350,7 @@ contains
scattxshg = cmfd%scattxs(h, g, i, j, k)
end if
! Negate the scattering xs
! Negate the scattering xs
val = -scattxshg
! Record value in matrix
@ -358,7 +358,7 @@ contains
end do SCATTR
end do ROWS
end do ROWS
! CSR requires n+1 row
call loss_matrix % new_row()

4
src/cmfd_prod_operator.F90
View file

@ -97,7 +97,7 @@ contains
end if
! loop around all other groups
! loop around all other groups
NFISS: do h = 1, ng
@ -122,7 +122,7 @@ contains
end do NFISS
end do ROWS
end do ROWS
! CSR requires n+1 row
call prod_matrix % new_row()

6
src/cmfd_solver.F90
View file

@ -65,7 +65,7 @@ contains
! Check for physical adjoint
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') &
physical_adjoint = .true.
physical_adjoint = .true.
! Start timer for build
call time_cmfdbuild % start()
@ -75,7 +75,7 @@ contains
! Check for mathematical adjoint calculation
if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') &
call compute_adjoint()
call compute_adjoint()
! Stop timer for build
call time_cmfdbuild % stop()
@ -286,7 +286,7 @@ contains
ROWS: do irow = 1, loss % n
COLS: do icol = loss % get_row(irow), loss % get_row(irow + 1) - 1
if (loss % get_col(icol) == prod % get_col(jcol) .and. &
jcol < prod % get_row(irow + 1)) then
jcol < prod % get_row(irow + 1)) then
loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol)
jcol = jcol + 1
end if

10
src/cross_section.F90
View file

@ -446,7 +446,7 @@ contains
! Calculate elastic cross section/factor
elastic = ZERO
if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then
elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, &
i_up)))
@ -455,7 +455,7 @@ contains
! Calculate fission cross section/factor
fission = ZERO
if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then
fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, &
i_up)))
@ -464,7 +464,7 @@ contains
! Calculate capture cross section/factor
capture = ZERO
if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. &
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then
capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, &
i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, &
i_up)))
@ -571,11 +571,11 @@ contains
! calculate interpolation factor
f = (E - nuc % energy_0K(i_grid)) &
& / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid))
& / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid))
! Calculate microscopic nuclide elastic cross section
xs_out = (ONE - f) * nuc % elastic_0K(i_grid) &
& + f * nuc % elastic_0K(i_grid + 1)
& + f * nuc % elastic_0K(i_grid + 1)
end function elastic_xs_0K

22
src/dict_header.F90
View file

@ -56,7 +56,7 @@ module dict_header
!===============================================================================
! 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.
! and DictIntInt has integer keys and values.
!===============================================================================
type, public :: DictCharInt
@ -105,7 +105,7 @@ contains
if (associated(elem)) then
elem % value = value
else
! Get hash
! Get hash
hash = dict_hash_key_ci(key)
! Create new element
@ -135,7 +135,7 @@ contains
if (associated(elem)) then
elem % value = value
else
! Get hash
! Get hash
hash = dict_hash_key_ii(key)
! Create new element
@ -156,13 +156,13 @@ contains
!===============================================================================
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))
@ -178,13 +178,13 @@ contains
end function dict_get_elem_ci
function dict_get_elem_ii(this, key) result(elem)
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))
@ -279,7 +279,7 @@ contains
character(*), intent(in) :: key
integer :: val
integer :: i
val = 0
@ -298,7 +298,7 @@ contains
integer, intent(in) :: key
integer :: val
val = 0
! Added the absolute val on val-1 since the sum in the do loop is
@ -420,7 +420,7 @@ contains
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
@ -440,5 +440,5 @@ contains
end if
end subroutine dict_clear_ii
end module dict_header

205
src/doppler.F90
View file

@ -52,62 +52,115 @@ contains
! Loop over incoming neutron energies
ENERGY_NEUTRON: do i = 1, n
sigma = ZERO
y = x(i)
y_sq = y*y
y_inv = ONE / y
y_inv_sq = y_inv / y
sigma = ZERO
y = x(i)
y_sq = y*y
y_inv = ONE / y
y_inv_sq = y_inv / y
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4
k = i
a = ZERO
call calculate_F(F_a, a)
k = i
a = ZERO
call calculate_F(F_a, a)
do while (a >= -4.0 .and. k > 1)
! Move to next point
F_b = F_a
k = k - 1
a = x(k) - y
do while (a >= -4.0 .and. k > 1)
! Move to next point
F_b = F_a
k = k - 1
a = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2)
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
! =======================================================================
! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE
if (k == 1 .and. a >= -4.0) then
! Since x = 0, this implies that a = -y
F_b = F_a
a = -y
if (k == 1 .and. a >= -4.0) then
! Since x = 0, this implies that a = -y
F_b = F_a
a = -y
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Calculate F and H functions
call calculate_F(F_a, a)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0))
end if
! Add contribution to broadened cross section
sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0))
end if
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
! =======================================================================
! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4
k = i
b = ZERO
call calculate_F(F_b, b)
k = i
b = ZERO
call calculate_F(F_b, b)
do while (b <= 4.0 .and. k < n)
do while (b <= 4.0 .and. k < n)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
! =======================================================================
! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE
if (k == n .and. b <= 4.0) then
! Calculate F function at last energy point
a = x(k) - y
call calculate_F(F_a, a)
! Add contribution to broadened cross section
sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0))
end if
! =======================================================================
! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4
if (y <= 4.0) then
! Swap signs on y
y = -y
y_inv = -y_inv
k = 1
! Calculate a and b based on 0 and x(1)
a = -y
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
call calculate_F(F_b, b)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0))
! Now progress forward doing the remainder of the second term
do while (b <= 4.0)
! Move to next point
F_a = F_b
k = k + 1
@ -119,69 +172,17 @@ contains
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) &
+ y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk
end do
sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk
end do
end if
! =======================================================================
! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE
if (k == n .and. b <= 4.0) then
! Calculate F function at last energy point
a = x(k) - y
call calculate_F(F_a, a)
! Add contribution to broadened cross section
sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0))
end if
! =======================================================================
! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4
if (y <= 4.0) then
! Swap signs on y
y = -y
y_inv = -y_inv
k = 1
! Calculate a and b based on 0 and x(1)
a = -y
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_a, a)
call calculate_F(F_b, b)
H = F_a - F_b
! Add contribution to broadened cross section
sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0))
! Now progress forward doing the remainder of the second term
do while (b <= 4.0)
! Move to next point
F_a = F_b
k = k + 1
b = x(k) - y
! Calculate F and H functions
call calculate_F(F_b, b)
H = F_a - F_b
! Calculate A(k), B(k), and slope terms
Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0)
Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0)
slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2)
! Add contribution to broadened cross section
sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk
end do
end if
! Set broadened cross section
sigmaNew(i) = sigma
! Set broadened cross section
sigmaNew(i) = sigma
end do ENERGY_NEUTRON

12
src/eigenvalue.F90
View file

@ -259,7 +259,7 @@ contains
! Write out source point if it's been specified for this batch
if ((sourcepoint_batch % contains(current_batch) .or. source_latest) .and. &
source_write) then
source_write) then
call write_source_point()
end if
@ -880,8 +880,8 @@ contains
!$omp do ordered schedule(static)
do i = 1, n_threads
!$omp ordered
master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
total = total + n_bank
master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank)
total = total + n_bank
!$omp end ordered
end do
!$omp end do
@ -891,10 +891,10 @@ contains
! Now copy the shared fission bank sites back to the master thread's copy.
if (thread_id == 0) then
n_bank = total
fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
n_bank = total
fission_bank(1:n_bank) = master_fission_bank(1:n_bank)
else
n_bank = 0
n_bank = 0
end if
!$omp end parallel

16
src/endf_header.F90
View file

@ -3,7 +3,7 @@ module endf_header
implicit none
!===============================================================================
! TAB1 represents a one-dimensional interpolable function
! TAB1 represents a one-dimensional interpolable function
!===============================================================================
type Tab1
@ -13,28 +13,28 @@ module endf_header
integer :: n_pairs ! # of pairs of (x,y) values
real(8), allocatable :: x(:) ! values of abscissa
real(8), allocatable :: y(:) ! values of ordinate
! Type-Bound procedures
contains
procedure :: clear => tab1_clear ! deallocates a Tab1 Object.
end type Tab1
contains
!===============================================================================
! TAB1_CLEAR deallocates the items in Tab1
!===============================================================================
subroutine tab1_clear(this)
class(Tab1), intent(inout) :: this ! The Tab1 to clear
if (allocated(this % nbt)) &
deallocate(this % nbt, this % int)
if (allocated(this % x)) &
deallocate(this % x, this % y)
end subroutine tab1_clear
end module endf_header

10
src/energy_grid.F90
View file

@ -92,15 +92,17 @@ contains
! Determine corresponding indices in nuclide grid to energies on
! equal-logarithmic grid
j = 1
do k = 0, M - 1
do k = 0, M
do while (log(nuc%energy(j + 1)/E_min) <= umesh(k))
! Ensure that for isotopes where maxval(nuc % energy) << E_max
! that there are no out-of-bounds issues.
if (j + 1 == nuc % n_grid) then
exit
end if
j = j + 1
end do
nuc % grid_index(k) = j
end do
! Set the last point explicitly so that we don't have out-of-bounds issues
nuc % grid_index(M) = size(nuc % energy) - 1
end do
deallocate(umesh)

69
src/geometry_header.F90
View file

@ -9,13 +9,13 @@ module geometry_header
!===============================================================================
type Universe
integer :: id ! Unique ID
integer :: type ! Type
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
integer :: id ! Unique ID
integer :: type ! Type
integer :: n_cells ! # of cells within
integer, allocatable :: cells(:) ! List of cells within
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
end type Universe
!===============================================================================
@ -119,14 +119,14 @@ module geometry_header
!===============================================================================
type Surface
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of surface
real(8), allocatable :: coeffs(:) ! Definition of surface
integer, allocatable :: &
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of surface
real(8), allocatable :: coeffs(:) ! Definition of surface
integer, allocatable :: &
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
end type Surface
!===============================================================================
@ -134,24 +134,29 @@ module geometry_header
!===============================================================================
type Cell
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of cell (normal, universe, lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: instances ! number of instances of this cell in the geom
integer :: material ! Material within cell (0 for universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: offset (:) ! Distribcell offset for tally counter
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
integer :: id ! Unique ID
character(len=52) :: name = "" ! User-defined name
integer :: type ! Type of cell (normal, universe,
! lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
integer :: instances ! number of instances of this cell in
! the geom
integer :: material ! Material within cell (0 for
! universe)
integer :: n_surfaces ! Number of surfaces within
integer, allocatable :: offset (:) ! Distribcell offset for tally
! counter
integer, allocatable :: &
& surfaces(:) ! List of surfaces bounding cell
! -- note that parentheses, union,
! etc operators will be listed here
! too
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)
real(8), allocatable :: rotation(:)
real(8), allocatable :: rotation_matrix(:,:)
! Rotation matrix and translation vector
real(8), allocatable :: translation(:)
real(8), allocatable :: rotation(:)
real(8), allocatable :: rotation_matrix(:,:)
end type Cell
! array index of universe 0

6
src/global.F90
View file

@ -529,9 +529,9 @@ contains
! Deallocate entropy mesh
if (associated(entropy_mesh)) then
if (allocated(entropy_mesh % lower_left)) &
deallocate(entropy_mesh % lower_left)
deallocate(entropy_mesh % lower_left)
if (allocated(entropy_mesh % upper_right)) &
deallocate(entropy_mesh % upper_right)
deallocate(entropy_mesh % upper_right)
if (allocated(entropy_mesh % width)) deallocate(entropy_mesh % width)
deallocate(entropy_mesh)
end if
@ -541,7 +541,7 @@ contains
if (associated(ufs_mesh)) then
if (allocated(ufs_mesh % lower_left)) deallocate(ufs_mesh % lower_left)
if (allocated(ufs_mesh % upper_right)) &
deallocate(ufs_mesh % upper_right)
deallocate(ufs_mesh % upper_right)
if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width)
deallocate(ufs_mesh)
end if

2
src/hdf5_interface.F90
View file

@ -7,7 +7,7 @@ module hdf5_interface
use, intrinsic :: ISO_C_BINDING
#ifdef MPI
use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL
use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL
#endif
implicit none

20
src/hdf5_summary.F90
View file

@ -186,7 +186,7 @@ contains
call su % write_data("lattice", "fill_type", &
group="geometry/cells/cell " // trim(to_str(c % id)))
call su % write_data(lattices(c % fill) % obj % id, "lattice", &
group="geometry/cells/cell " // trim(to_str(c % id)))
group="geometry/cells/cell " // trim(to_str(c % id)))
end select
! Write list of bounding surfaces
@ -373,7 +373,7 @@ contains
! Write lattice universes.
allocate(lattice_universes(lat % n_cells(1), lat % n_cells(2), &
&lat % n_cells(3)))
&lat % n_cells(3)))
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
@ -399,13 +399,13 @@ contains
! Write lattice center, pitch and outer universe.
if (lat % is_3d) then
call su % write_data(lat % center, "center", length=3, &
call su % write_data(lat % center, "center", length=3, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
else
call su % write_data(lat % center, "center", length=2, &
call su % write_data(lat % center, "center", length=2, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
end if
if (lat % is_3d) then
call su % write_data(lat % pitch, "pitch", length=2, &
group="geometry/lattices/lattice " // trim(to_str(lat % id)))
@ -429,7 +429,7 @@ contains
! Write lattice universes.
allocate(lattice_universes(2*lat % n_rings - 1, 2*lat % n_rings - 1, &
&lat % n_axial))
&lat % n_axial))
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
@ -647,12 +647,12 @@ contains
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYIN)
call su % write_data("energy", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
case(FILTER_ENERGYOUT)
call su % write_data("energyout", "type_name", &
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
group="tallies/tally " // trim(to_str(t % id)) &
// "/filter " // trim(to_str(j)))
end select
end do FILTER_LOOP

4
src/initialize.F90
View file

@ -916,9 +916,9 @@ contains
thread_id = omp_get_thread_num()
if (thread_id == 0) then
allocate(fission_bank(3*work))
allocate(fission_bank(3*work))
else
allocate(fission_bank(3*work/n_threads))
allocate(fission_bank(3*work/n_threads))
end if
!$omp end parallel
allocate(master_fission_bank(3*work), STAT=alloc_err)

238
src/input_xml.F90
View file

@ -472,7 +472,7 @@ contains
! Check for type of energy distribution
type = ''
if (check_for_node(node_dist, "type")) &
call get_node_value(node_dist, "type", type)
call get_node_value(node_dist, "type", type)
select case (to_lower(type))
case ('monoenergetic')
external_source % type_energy = SRC_ENERGY_MONO
@ -798,7 +798,7 @@ contains
if (.not. source_separate) then
do i = 1, n_source_points
if (.not. statepoint_batch % contains(sourcepoint_batch % &
get_item(i))) then
get_item(i))) then
call fatal_error('Sourcepoint batches are not a subset&
& of statepoint batches.')
end if
@ -811,7 +811,7 @@ contains
call get_node_value(doc, "no_reduce", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
reduce_tallies = .false.
reduce_tallies = .false.
end if
! Check if the user has specified to use confidence intervals for
@ -884,11 +884,11 @@ contains
&// trim(to_str(i)) // " in settings.xml file!")
end if
call get_node_value(node_scatterer, "nuclide", &
nuclides_0K(i) % nuclide)
nuclides_0K(i) % nuclide)
if (check_for_node(node_scatterer, "method")) then
call get_node_value(node_scatterer, "method", &
nuclides_0K(i) % scheme)
nuclides_0K(i) % scheme)
end if
! check to make sure xs name for which method is applied is given
@ -898,7 +898,7 @@ contains
&// " given in cross_sections.xml")
end if
call get_node_value(node_scatterer, "xs_label", &
nuclides_0K(i) % name)
nuclides_0K(i) % name)
! check to make sure 0K xs name for which method is applied is given
if (.not. check_for_node(node_scatterer, "xs_label_0K")) then
@ -906,11 +906,11 @@ contains
&// trim(to_str(i)) // " given in cross_sections.xml")
end if
call get_node_value(node_scatterer, "xs_label_0K", &
nuclides_0K(i) % name_0K)
nuclides_0K(i) % name_0K)
if (check_for_node(node_scatterer, "E_min")) then
call get_node_value(node_scatterer, "E_min", &
nuclides_0K(i) % E_min)
nuclides_0K(i) % E_min)
end if
! check that E_min is non-negative
@ -921,7 +921,7 @@ contains
if (check_for_node(node_scatterer, "E_max")) then
call get_node_value(node_scatterer, "E_max", &
nuclides_0K(i) % E_max)
nuclides_0K(i) % E_max)
end if
! check that E_max is not less than E_min
@ -1081,7 +1081,7 @@ contains
! Read material
word = ''
if (check_for_node(node_cell, "material")) &
call get_node_value(node_cell, "material", word)
call get_node_value(node_cell, "material", word)
select case(to_lower(word))
case ('void')
c % material = MATERIAL_VOID
@ -1248,7 +1248,7 @@ contains
! Copy and interpret surface type
word = ''
if (check_for_node(node_surf, "type")) &
call get_node_value(node_surf, "type", word)
call get_node_value(node_surf, "type", word)
select case(to_lower(word))
case ('x-plane')
s % type = SURF_PX
@ -1306,7 +1306,7 @@ contains
! Boundary conditions
word = ''
if (check_for_node(node_surf, "boundary")) &
call get_node_value(node_surf, "boundary", word)
call get_node_value(node_surf, "boundary", word)
select case (to_lower(word))
case ('transmission', 'transmit', '')
s % bc = BC_TRANSMIT
@ -1447,7 +1447,7 @@ contains
do k = 0, n_y - 1
do j = 1, n_x
lat % universes(j, n_y - k, m) = &
&temp_int_array(j + n_x*k + n_x*n_y*(m-1))
&temp_int_array(j + n_x*k + n_x*n_y*(m-1))
end do
end do
end do
@ -1713,7 +1713,7 @@ contains
! Copy default cross section if present
if (check_for_node(doc, "default_xs")) &
call get_node_value(doc, "default_xs", default_xs)
call get_node_value(doc, "default_xs", default_xs)
! Get pointer to list of XML <material>
call get_node_list(doc, "material", node_mat_list)
@ -1844,7 +1844,7 @@ contains
! store full name
call get_node_value(node_nuc, "name", temp_str)
if (check_for_node(node_nuc, "xs")) &
call get_node_value(node_nuc, "xs", name)
call get_node_value(node_nuc, "xs", name)
name = trim(temp_str) // "." // trim(name)
! save name and density to list
@ -1853,7 +1853,7 @@ contains
! Check if no atom/weight percents were specified or if both atom and
! weight percents were specified
if (.not.check_for_node(node_nuc, "ao") .and. &
.not.check_for_node(node_nuc, "wo")) then
.not.check_for_node(node_nuc, "wo")) then
call fatal_error("No atom or weight percent specified for nuclide " &
&// trim(name))
elseif (check_for_node(node_nuc, "ao") .and. &
@ -1903,7 +1903,7 @@ contains
! Check if no atom/weight percents were specified or if both atom and
! weight percents were specified
if (.not.check_for_node(node_ele, "ao") .and. &
.not.check_for_node(node_ele, "wo")) then
.not.check_for_node(node_ele, "wo")) then
call fatal_error("No atom or weight percent specified for element " &
&// trim(name))
elseif (check_for_node(node_ele, "ao") .and. &
@ -2010,7 +2010,7 @@ contains
! Determine name of S(a,b) table
if (.not.check_for_node(node_sab, "name") .or. &
.not.check_for_node(node_sab, "xs")) then
.not.check_for_node(node_sab, "xs")) then
call fatal_error("Need to specify <name> and <xs> for S(a,b) &
&table.")
end if
@ -2157,7 +2157,7 @@ contains
call get_node_value(doc, "assume_separate", temp_str)
temp_str = to_lower(temp_str)
if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') &
assume_separate = .true.
assume_separate = .true.
end if
! ==========================================================================
@ -2185,7 +2185,7 @@ contains
! Read mesh type
temp_str = ''
if (check_for_node(node_mesh, "type")) &
call get_node_value(node_mesh, "type", temp_str)
call get_node_value(node_mesh, "type", temp_str)
select case (to_lower(temp_str))
case ('rect', 'rectangle', 'rectangular')
m % type = LATTICE_RECT
@ -2227,14 +2227,14 @@ contains
! Make sure both upper-right or width were specified
if (check_for_node(node_mesh, "upper_right") .and. &
check_for_node(node_mesh, "width")) then
check_for_node(node_mesh, "width")) then
call fatal_error("Cannot specify both <upper_right> and <width> on a &
&tally mesh.")
end if
! Make sure either upper-right or width was specified
if (.not.check_for_node(node_mesh, "upper_right") .and. &
.not.check_for_node(node_mesh, "width")) then
.not.check_for_node(node_mesh, "width")) then
call fatal_error("Must specify either <upper_right> and <width> on a &
&tally mesh.")
end if
@ -2242,7 +2242,7 @@ contains
if (check_for_node(node_mesh, "width")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "width") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <width> must be the same as &
&the number of entries on <lower_left>.")
end if
@ -2260,7 +2260,7 @@ contains
elseif (check_for_node(node_mesh, "upper_right")) then
! Check to ensure width has same dimensions
if (get_arraysize_double(node_mesh, "upper_right") /= &
get_arraysize_double(node_mesh, "lower_left")) then
get_arraysize_double(node_mesh, "lower_left")) then
call fatal_error("Number of entries on <upper_right> must be the &
&same as the number of entries on <lower_left>.")
end if
@ -2323,7 +2323,7 @@ contains
! Copy tally name
if (check_for_node(node_tal, "name")) &
call get_node_value(node_tal, "name", t % name)
call get_node_value(node_tal, "name", t % name)
! =======================================================================
! READ DATA FOR FILTERS
@ -2345,13 +2345,13 @@ contains
! Convert filter type to lower case
temp_str = ''
if (check_for_node(node_filt, "type")) &
call get_node_value(node_filt, "type", temp_str)
call get_node_value(node_filt, "type", temp_str)
temp_str = to_lower(temp_str)
! Determine number of bins
if (check_for_node(node_filt, "bins")) then
if (trim(temp_str) == 'energy' .or. &
trim(temp_str) == 'energyout') then
trim(temp_str) == 'energyout') then
n_words = get_arraysize_double(node_filt, "bins")
else
n_words = get_arraysize_integer(node_filt, "bins")
@ -2618,7 +2618,7 @@ contains
if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
@ -2661,7 +2661,7 @@ contains
if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
@ -2684,7 +2684,7 @@ contains
if (starts_with(score_name,trim(MOMENT_N_STRS(imomstr)))) then
n_order_pos = scan(score_name,'0123456789')
n_order = int(str_to_int( &
score_name(n_order_pos:(len_trim(score_name)))),4)
score_name(n_order_pos:(len_trim(score_name)))),4)
if (n_order > MAX_ANG_ORDER) then
! User requested too many orders; throw a warning and set to the
! maximum order.
@ -3006,7 +3006,7 @@ contains
! Get the trigger type - "variance", "std_dev" or "rel_err"
if (check_for_node(node_trigger, "type")) then
call get_node_value(node_trigger, "type", temp_str)
temp_str = to_lower(temp_str)
temp_str = to_lower(temp_str)
else
call fatal_error("Must specify trigger type for tally " // &
trim(to_str(t % id)) // " in tally XML file.")
@ -3105,7 +3105,7 @@ contains
! Increment the overall trigger index
trig_ind = trig_ind + 1
end if
end if
end do SCORE_LOOP
! Deallocate the list of tally scores used to create triggers
@ -3223,7 +3223,7 @@ contains
! Copy plot type
temp_str = 'slice'
if (check_for_node(node_plot, "type")) &
call get_node_value(node_plot, "type", temp_str)
call get_node_value(node_plot, "type", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("slice")
@ -3238,7 +3238,7 @@ contains
! Set output file path
filename = trim(to_str(pl % id)) // "_plot"
if (check_for_node(node_plot, "filename")) &
call get_node_value(node_plot, "filename", filename)
call get_node_value(node_plot, "filename", filename)
select case (pl % type)
case (PLOT_TYPE_SLICE)
pl % path_plot = trim(path_input) // trim(filename) // ".ppm"
@ -3283,7 +3283,7 @@ contains
if (pl % type == PLOT_TYPE_SLICE) then
temp_str = 'xy'
if (check_for_node(node_plot, "basis")) &
call get_node_value(node_plot, "basis", temp_str)
call get_node_value(node_plot, "basis", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("xy")
@ -3338,7 +3338,7 @@ contains
! Copy plot color type and initialize all colors randomly
temp_str = "cell"
if (check_for_node(node_plot, "color")) &
call get_node_value(node_plot, "color", temp_str)
call get_node_value(node_plot, "color", temp_str)
temp_str = to_lower(temp_str)
select case (trim(temp_str))
case ("cell")
@ -3452,7 +3452,7 @@ contains
! Ensure that there is a linewidth for this meshlines specification
if (check_for_node(node_meshlines, "linewidth")) then
call get_node_value(node_meshlines, "linewidth", &
pl % meshlines_width)
pl % meshlines_width)
else
call fatal_error("Must specify a linewidth for meshlines &
&specification in plot " // trim(to_str(pl % id)))
@ -3468,7 +3468,7 @@ contains
end if
call get_node_array(node_meshlines, "color", &
pl % meshlines_color % rgb)
pl % meshlines_color % rgb)
else
pl % meshlines_color % rgb = (/ 0, 0, 0 /)
@ -3494,8 +3494,8 @@ contains
end if
i_mesh = cmfd_tallies(1) % &
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % &
int_bins(1)
filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % &
int_bins(1)
pl % meshlines_mesh => meshes(i_mesh)
case ('entropy')
@ -3654,9 +3654,9 @@ contains
! Check if cross_sections.xml exists
inquire(FILE=path_cross_sections, EXIST=file_exists)
if (.not. file_exists) then
! Could not find cross_sections.xml file
call fatal_error("Cross sections XML file '" &
&// trim(path_cross_sections) // "' does not exist!")
! Could not find cross_sections.xml file
call fatal_error("Cross sections XML file '" &
&// trim(path_cross_sections) // "' does not exist!")
end if
call write_message("Reading cross sections XML file...", 5)
@ -3665,28 +3665,28 @@ contains
call open_xmldoc(doc, path_cross_sections)
if (check_for_node(doc, "directory")) then
! Copy directory information if present
call get_node_value(doc, "directory", directory)
! Copy directory information if present
call get_node_value(doc, "directory", directory)
else
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
! If no directory is listed in cross_sections.xml, by default select the
! directory in which the cross_sections.xml file resides
i = index(path_cross_sections, "/", BACK=.true.)
directory = path_cross_sections(1:i)
end if
! determine whether binary/ascii
temp_str = ''
if (check_for_node(doc, "filetype")) &
call get_node_value(doc, "filetype", temp_str)
call get_node_value(doc, "filetype", temp_str)
if (trim(temp_str) == 'ascii') then
filetype = ASCII
filetype = ASCII
elseif (trim(temp_str) == 'binary') then
filetype = BINARY
filetype = BINARY
elseif (len_trim(temp_str) == 0) then
filetype = ASCII
filetype = ASCII
else
call fatal_error("Unknown filetype in cross_sections.xml: " &
&// trim(temp_str))
call fatal_error("Unknown filetype in cross_sections.xml: " &
&// trim(temp_str))
end if
! copy default record length and entries for binary files
@ -3701,83 +3701,83 @@ contains
! Allocate xs_listings array
if (n_listings == 0) then
call fatal_error("No ACE table listings present in cross_sections.xml &
&file!")
call fatal_error("No ACE table listings present in cross_sections.xml &
&file!")
else
allocate(xs_listings(n_listings))
allocate(xs_listings(n_listings))
end if
do i = 1, n_listings
listing => xs_listings(i)
listing => xs_listings(i)
! Get pointer to ace table XML node
call get_list_item(node_ace_list, i, node_ace)
! Get pointer to ace table XML node
call get_list_item(node_ace_list, i, node_ace)
! copy a number of attributes
call get_node_value(node_ace, "name", listing % name)
if (check_for_node(node_ace, "alias")) &
call get_node_value(node_ace, "alias", listing % alias)
call get_node_value(node_ace, "zaid", listing % zaid)
call get_node_value(node_ace, "awr", listing % awr)
if (check_for_node(node_ace, "temperature")) &
call get_node_value(node_ace, "temperature", listing % kT)
call get_node_value(node_ace, "location", listing % location)
! copy a number of attributes
call get_node_value(node_ace, "name", listing % name)
if (check_for_node(node_ace, "alias")) &
call get_node_value(node_ace, "alias", listing % alias)
call get_node_value(node_ace, "zaid", listing % zaid)
call get_node_value(node_ace, "awr", listing % awr)
if (check_for_node(node_ace, "temperature")) &
call get_node_value(node_ace, "temperature", listing % kT)
call get_node_value(node_ace, "location", listing % location)
! determine type of cross section
if (ends_with(listing % name, 'c')) then
listing % type = ACE_NEUTRON
elseif (ends_with(listing % name, 't')) then
listing % type = ACE_THERMAL
end if
! determine type of cross section
if (ends_with(listing % name, 'c')) then
listing % type = ACE_NEUTRON
elseif (ends_with(listing % name, 't')) then
listing % type = ACE_THERMAL
end if
! set filetype, record length, and number of entries
if (check_for_node(node_ace, "filetype")) then
temp_str = ''
call get_node_value(node_ace, "filetype", temp_str)
if (temp_str == 'ascii') then
listing % filetype = ASCII
else if (temp_str == 'binary') then
listing % filetype = BINARY
end if
else
listing % filetype = filetype
end if
! set filetype, record length, and number of entries
if (check_for_node(node_ace, "filetype")) then
temp_str = ''
call get_node_value(node_ace, "filetype", temp_str)
if (temp_str == 'ascii') then
listing % filetype = ASCII
else if (temp_str == 'binary') then
listing % filetype = BINARY
end if
else
listing % filetype = filetype
end if
! Set record length and entries for binary files
if (filetype == BINARY) then
listing % recl = recl
listing % entries = entries
end if
! Set record length and entries for binary files
if (filetype == BINARY) then
listing % recl = recl
listing % entries = entries
end if
! determine metastable state
if (.not.check_for_node(node_ace, "metastable")) then
listing % metastable = .false.
else
listing % metastable = .true.
end if
! determine metastable state
if (.not.check_for_node(node_ace, "metastable")) then
listing % metastable = .false.
else
listing % metastable = .true.
end if
! determine path of cross section table
if (check_for_node(node_ace, "path")) then
call get_node_value(node_ace, "path", temp_str)
else
call fatal_error("Path missing for isotope " // listing % name)
end if
! determine path of cross section table
if (check_for_node(node_ace, "path")) then
call get_node_value(node_ace, "path", temp_str)
else
call fatal_error("Path missing for isotope " // listing % name)
end if
if (starts_with(temp_str, '/')) then
listing % path = trim(temp_str)
else
if (ends_with(directory,'/')) then
listing % path = trim(directory) // trim(temp_str)
else
listing % path = trim(directory) // '/' // trim(temp_str)
end if
end if
if (starts_with(temp_str, '/')) then
listing % path = trim(temp_str)
else
if (ends_with(directory,'/')) then
listing % path = trim(directory) // trim(temp_str)
else
listing % path = trim(directory) // '/' // trim(temp_str)
end if
end if
! create dictionary entry for both name and alias
call xs_listing_dict % add_key(to_lower(listing % name), i)
if (check_for_node(node_ace, "alias")) then
call xs_listing_dict % add_key(to_lower(listing % alias), i)
end if
! create dictionary entry for both name and alias
call xs_listing_dict % add_key(to_lower(listing % name), i)
if (check_for_node(node_ace, "alias")) then
call xs_listing_dict % add_key(to_lower(listing % alias), i)
end if
end do
! Check that 0K nuclides are listed in the cross_sections.xml file

18
src/list_header.F90
View file

@ -45,7 +45,7 @@ module list_header
private
integer :: count = 0 ! Number of elements in list
! Used in get_item for fast sequential lookups
! Used in get_item for fast sequential lookups
integer :: last_index = huge(0)
type(ListElemInt), pointer :: last_elem => null()
@ -67,7 +67,7 @@ module list_header
private
integer :: count = 0 ! Number of elements in list
! Used in get_item for fast sequential lookups
! Used in get_item for fast sequential lookups
integer :: last_index = huge(0)
type(ListElemReal), pointer :: last_elem => null()
@ -89,7 +89,7 @@ module list_header
private
integer :: count = 0 ! Number of elements in list
! Used in get_item for fast sequential lookups
! Used in get_item for fast sequential lookups
integer :: last_index = huge(0)
type(ListElemChar), pointer :: last_elem => null()
@ -198,7 +198,7 @@ contains
type(ListElemInt), pointer :: current => null()
type(ListElemInt), pointer :: next => null()
if (this % count > 0) then
current => this % head
do while (associated(current))
@ -224,7 +224,7 @@ contains
type(ListElemReal), pointer :: current => null()
type(ListElemReal), pointer :: next => null()
if (this % count > 0) then
current => this % head
do while (associated(current))
@ -250,7 +250,7 @@ contains
type(ListElemChar), pointer :: current => null()
type(ListElemChar), pointer :: next => null()
if (this % count > 0) then
current => this % head
do while (associated(current))
@ -518,7 +518,7 @@ contains
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % prev => elem % prev
new_elem % next => elem
@ -573,7 +573,7 @@ contains
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % prev => elem % prev
new_elem % next => elem
@ -628,7 +628,7 @@ contains
! Allocate new element
allocate(new_elem)
new_elem % data = data
! Put it before the i-th element
new_elem % prev => elem % prev
new_elem % next => elem

241
src/math.F90
View file

@ -143,20 +143,20 @@ contains
pnx = 7.875_8 * (x ** 5) - 8.75_8 * x * x * x + 1.875 * x
case(6)
pnx = 14.4375_8 * (x ** 6) - 19.6875_8 * (x ** 4) + &
6.5625_8 * x * x - 0.3125_8
6.5625_8 * x * x - 0.3125_8
case(7)
pnx = 26.8125_8 * (x ** 7) - 43.3125_8 * (x ** 5) + &
19.6875_8 * x * x * x - 2.1875_8 * x
19.6875_8 * x * x * x - 2.1875_8 * x
case(8)
pnx = 50.2734375_8 * (x ** 8) - 93.84375_8 * (x ** 6) + &
54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8
54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8
case(9)
pnx = 94.9609375_8 * (x ** 9) - 201.09375_8 * (x ** 7) + &
140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x
140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x
case(10)
pnx = 180.42578125_8 * (x ** 10) - 427.32421875_8 * (x ** 8) + &
351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + &
13.53515625_8 * x * x - 0.24609375_8
351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + &
13.53515625_8 * x * x - 0.24609375_8
case default
pnx = ONE ! correct for case(0), incorrect for the rest
end select
@ -215,12 +215,12 @@ contains
rn(2) = 1.93649167310371_8 * w*(w2m1) * sin(TWO*phi)
! l = 3, m = -1
rn(3) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
sin(phi)
sin(phi)
! l = 3, m = 0
rn(4) = 2.5_8 * w**3 - 1.5_8 * w
! l = 3, m = 1
rn(5) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * &
cos(phi)
cos(phi)
! l = 3, m = 2
rn(6) = 1.93649167310371_8 * w*(w2m1) * cos(TWO*phi)
! l = 3, m = 3
@ -232,18 +232,18 @@ contains
rn(2) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * sin(THREE* phi)
! l = 4, m = -2
rn(3) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
sin(TWO*phi)
sin(TWO*phi)
! l = 4, m = -1
rn(4) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * &
sin(phi)
rn(4) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)&
* sin(phi)
! l = 4, m = 0
rn(5) = 4.375_8 * w**4 - 3.75_8 * w**2 + 0.375_8
! l = 4, m = 1
rn(6) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * &
cos(phi)
rn(6) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)&
* cos(phi)
! l = 4, m = 2
rn(7) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * &
cos(TWO*phi)
cos(TWO*phi)
! l = 4, m = 3
rn(8) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * cos(THREE* phi)
! l = 4, m = 4
@ -255,24 +255,26 @@ contains
rn(2) = 2.21852991866236_8 * w*(w2m1)**2 * sin(4.0_8*phi)
! l = 5, m = -3
rn(3) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi)
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi)
! l = 5, m = -2
rn(4) = 0.0487950036474267_8 * (w2m1)*((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * &
sin(TWO*phi)
rn(4) = 0.0487950036474267_8 * (w2m1) &
* ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * sin(TWO*phi)
! l = 5, m = -1
rn(5) = 0.258198889747161_8*sqrt(w2m1)* &
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * sin(phi)
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) &
* sin(phi)
! l = 5, m = 0
rn(6) = 7.875_8 * w**5 - 8.75_8 * w**3 + 1.875_8 * w
! l = 5, m = 1
rn(7) = 0.258198889747161_8*sqrt(w2m1)* &
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * cos(phi)
((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) &
* cos(phi)
! l = 5, m = 2
rn(8) = 0.0487950036474267_8 * (w2m1)* &
((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi)
((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi)
! l = 5, m = 3
rn(9) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* &
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi)
((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi)
! l = 5, m = 4
rn(10) = 2.21852991866236_8 * w*(w2m1)**2 * cos(4.0_8*phi)
! l = 5, m = 5
@ -284,30 +286,34 @@ contains
rn(2) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * sin(5.0_8*phi)
! l = 6, m = -4
rn(3) = 0.00104990131391452_8 * (w2m1)**2 * &
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi)
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi)
! l = 6, m = -3
rn(4) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi)
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi)
! l = 6, m = -2
rn(5) = 0.0345032779671177_8 * (w2m1) * &
((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * sin(TWO*phi)
((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) &
* sin(TWO*phi)
! l = 6, m = -1
rn(6) = 0.218217890235992_8*sqrt(w2m1) * &
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * sin(phi)
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) &
* sin(phi)
! l = 6, m = 0
rn(7) = 14.4375_8 * w**6 - 19.6875_8 * w**4 + 6.5625_8 * w**2 - 0.3125_8
! l = 6, m = 1
rn(8) = 0.218217890235992_8*sqrt(w2m1) * &
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * cos(phi)
((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) &
* cos(phi)
! l = 6, m = 2
rn(9) = 0.0345032779671177_8 * (w2m1) * &
((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * cos(TWO*phi)
((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) &
* cos(TWO*phi)
! l = 6, m = 3
rn(10) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * &
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi)
((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi)
! l = 6, m = 4
rn(11) = 0.00104990131391452_8 * (w2m1)**2 * &
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi)
((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi)
! l = 6, m = 5
rn(12) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * cos(5.0_8*phi)
! l = 6, m = 6
@ -319,42 +325,43 @@ contains
rn(2) = 2.42182459624969_8 * w*(w2m1)**3 * sin(6.0_8*phi)
! l = 7, m = -5
rn(3) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* &
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi)
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi)
! l = 7, m = -4
rn(4) = 0.000548293079133141_8 * (w2m1)**2* &
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi)
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi)
! l = 7, m = -3
rn(5) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
sin(THREE*phi)
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
sin(THREE*phi)
! l = 7, m = -2
rn(6) = 0.025717224993682_8 * (w2m1)* &
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
sin(TWO*phi)
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
sin(TWO*phi)
! l = 7, m = -1
rn(7) = 0.188982236504614_8*sqrt(w2m1)* &
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi)
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi)
! l = 7, m = 0
rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 * w
rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 &
* w
! l = 7, m = 1
rn(9) = 0.188982236504614_8*sqrt(w2m1)* &
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi)
((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + &
(945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi)
! l = 7, m = 2
rn(10) = 0.025717224993682_8 * (w2m1)* &
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
cos(TWO*phi)
((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* &
cos(TWO*phi)
! l = 7, m = 3
rn(11) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* &
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
cos(THREE*phi)
((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* &
cos(THREE*phi)
! l = 7, m = 4
rn(12) = 0.000548293079133141_8 * (w2m1)**2 * &
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi)
((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi)
! l = 7, m = 5
rn(13) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* &
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi)
((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi)
! l = 7, m = 6
rn(14) = 2.42182459624969_8 * w*(w2m1)**3 * cos(6.0_8*phi)
! l = 7, m = 7
@ -366,48 +373,50 @@ contains
rn(2) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * sin(7.0_8*phi)
! l = 8, m = -6
rn(3) = 6.77369783729086d-6*(w2m1)**3* &
((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi)
((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi)
! l = 8, m = -5
rn(4) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)* &
((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi)
((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi)
! l = 8, m = -4
rn(5) = 0.000316557156832328_8 * (w2m1)**2* &
((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * sin(4.0_8*phi)
((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 &
+ 10395.0_8/8.0_8) * sin(4.0_8*phi)
! l = 8, m = -3
rn(6) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + (10395.0_8/8.0_8)*w) * sin(THREE*phi)
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 &
+ (10395.0_8/8.0_8)*w) * sin(THREE*phi)
! l = 8, m = -2
rn(7) = 0.0199204768222399_8 * (w2m1)* &
((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + &
(10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi)
((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + &
(10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi)
! l = 8, m = -1
rn(8) = 0.166666666666667_8*sqrt(w2m1)* &
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi)
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi)
! l = 8, m = 0
rn(9) = 50.2734375_8 * w**8 - 93.84375_8 * w**6 + 54.140625_8 * w**4 -&
9.84375_8 * w**2 + 0.2734375_8
9.84375_8 * w**2 + 0.2734375_8
! l = 8, m = 1
rn(10) = 0.166666666666667_8*sqrt(w2m1)* &
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi)
((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + &
(3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi)
! l = 8, m = 2
rn(11) = 0.0199204768222399_8 * (w2m1)*((45045.0_8/16.0_8)*w**6- &
45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - &
315.0_8/16.0_8) * cos(TWO*phi)
45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - &
315.0_8/16.0_8) * cos(TWO*phi)
! l = 8, m = 3
rn(12) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* &
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + &
(10395.0_8/8.0_8)*w) * cos(THREE*phi)
((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + &
(10395.0_8/8.0_8)*w) * cos(THREE*phi)
! l = 8, m = 4
rn(13) = 0.000316557156832328_8 * (w2m1)**2*((675675.0_8/8.0_8)*w**4 - &
135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi)
135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi)
! l = 8, m = 5
rn(14) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 - &
135135.0_8/TWO*w) * cos(5.0_8*phi)
rn(14) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 -&
135135.0_8/TWO*w) * cos(5.0_8*phi)
! l = 8, m = 6
rn(15) = 6.77369783729086d-6*(w2m1)**3*((2027025.0_8/TWO)*w**2 - &
135135.0_8/TWO) * cos(6.0_8*phi)
135135.0_8/TWO) * cos(6.0_8*phi)
! l = 8, m = 7
rn(16) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * cos(7.0_8*phi)
! l = 8, m = 8
@ -419,54 +428,56 @@ contains
rn(2) = 2.58397773170915_8 * w*(w2m1)**4 * sin(8.0_8*phi)
! l = 9, m = -7
rn(3) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* &
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi)
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi)
! l = 9, m = -6
rn(4) = 3.02928976464514d-6*(w2m1)**3* &
((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi)
((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi)
! l = 9, m = -5
rn(5) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)* &
((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + &
135135.0_8/8.0_8) * sin(5.0_8*phi)
((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + &
135135.0_8/8.0_8) * sin(5.0_8*phi)
! l = 9, m = -4
rn(6) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi)
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi)
! l = 9, m = -3
rn(7) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)* &
((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + &
(135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi)
((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + &
(135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi)
! l = 9, m = -2
rn(8) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7- &
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/16.0_8 * w)* &
sin(TWO*phi)
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 &
- 3465.0_8/16.0_8 * w) * sin(TWO*phi)
! l = 9, m = -1
rn(9) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - &
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * sin(phi)
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 &
* w**2 + 315.0_8/128.0_8) * sin(phi)
! l = 9, m = 0
rn(10) = 94.9609375_8 * w**9 - 201.09375_8 * w**7 + 140.765625_8 * w**5- &
36.09375_8 * w**3 + 2.4609375_8 * w
36.09375_8 * w**3 + 2.4609375_8 * w
! l = 9, m = 1
rn(11) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - &
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * cos(phi)
45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 &
* w**2 + 315.0_8/128.0_8) * cos(phi)
! l = 9, m = 2
rn(12) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7 - &
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/ 16.0_8 * w) * &
cos(TWO*phi)
135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 &
- 3465.0_8/ 16.0_8 * w) * cos(TWO*phi)
! l = 9, m = 3
rn(13) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)*w**6 - &
675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8)* &
cos(THREE*phi)
rn(13) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)&
*w**6 - 675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 &
- 3465.0_8/16.0_8)* cos(THREE*phi)
! l = 9, m = 4
rn(14) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - &
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi)
675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi)
! l = 9, m = 5
rn(15) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)*((11486475.0_8/8.0_8)* &
w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi)
w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi)
! l = 9, m = 6
rn(16) = 3.02928976464514d-6*(w2m1)**3*((11486475.0_8/TWO)*w**3 - &
2027025.0_8/TWO*w) * cos(6.0_8*phi)
2027025.0_8/TWO*w) * cos(6.0_8*phi)
! l = 9, m = 7
rn(17) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* &
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi)
((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi)
! l = 9, m = 8
rn(18) = 2.58397773170915_8 * w*(w2m1)**4 * cos(8.0_8*phi)
! l = 9, m = 9
@ -478,65 +489,65 @@ contains
rn(2) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * sin(9.0_8*phi)
! l = 10, m = -8
rn(3) = 2.49953651452314d-8*(w2m1)**4*((654729075.0_8/TWO)*w**2 - &
34459425.0_8/TWO) * sin(8.0_8*phi)
34459425.0_8/TWO) * sin(8.0_8*phi)
! l = 10, m = -7
rn(4) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* &
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi)
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi)
! l = 10, m = -6
rn(5) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - &
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi)
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi)
! l = 10, m = -5
rn(6) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* &
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * sin(5.0_8*phi)
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * sin(5.0_8*phi)
! l = 10, m = -4
rn(7) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - &
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * sin(4.0_8*phi)
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * sin(4.0_8*phi)
! l = 10, m = -3
rn(8) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi)
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi)
! l = 10, m = -2
rn(9) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - &
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - &
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi)
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - &
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi)
! l = 10, m = -1
rn(10) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - &
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - &
15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi)
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - &
15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi)
! l = 10, m = 0
rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 * w**6 - &
117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8
rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 &
* w**6 - 117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8
! l = 10, m = 1
rn(12) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - &
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ &
32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi)
109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ &
32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi)
! l = 10, m = 2
rn(13) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - &
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -&
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi)
765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -&
45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi)
! l = 10, m = 3
rn(14) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* &
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi)
((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + &
(675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi)
! l = 10, m = 4
rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - &
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * cos(4.0_8*phi)
rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 -&
11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - &
45045.0_8/16.0_8) * cos(4.0_8*phi)
! l = 10, m = 5
rn(16) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* &
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * cos(5.0_8*phi)
((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + &
(2027025.0_8/8.0_8)*w) * cos(5.0_8*phi)
! l = 10, m = 6
rn(17) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - &
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi)
34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi)
! l = 10, m = 7
rn(18) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* &
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi)
((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi)
! l = 10, m = 8
rn(19) = 2.49953651452314d-8*(w2m1)**4* &
((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi)
((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi)
! l = 10, m = 9
rn(20) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * cos(9.0_8*phi)
! l = 10, m = 10

26
src/matrix_header.F90
View file

@ -11,19 +11,19 @@ module matrix_header
integer, allocatable :: row(:) ! csr row vector
integer, allocatable :: col(:) ! column vector
real(8), allocatable :: val(:) ! matrix value vector
contains
procedure :: create => matrix_create
procedure :: destroy => matrix_destroy
procedure :: add_value => matrix_add_value
procedure :: new_row => matrix_new_row
procedure :: assemble => matrix_assemble
procedure :: get_row => matrix_get_row
procedure :: get_col => matrix_get_col
procedure :: vector_multiply => matrix_vector_multiply
procedure :: search_indices => matrix_search_indices
procedure :: write => matrix_write
procedure :: copy => matrix_copy
procedure :: transpose => matrix_transpose
contains
procedure :: create => matrix_create
procedure :: destroy => matrix_destroy
procedure :: add_value => matrix_add_value
procedure :: new_row => matrix_new_row
procedure :: assemble => matrix_assemble
procedure :: get_row => matrix_get_row
procedure :: get_col => matrix_get_col
procedure :: vector_multiply => matrix_vector_multiply
procedure :: search_indices => matrix_search_indices
procedure :: write => matrix_write
procedure :: copy => matrix_copy
procedure :: transpose => matrix_transpose
end type matrix
contains

28
src/mpiio_interface.F90
View file

@ -129,13 +129,13 @@ contains
integer, intent(inout) :: buffer ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_integer
@ -341,13 +341,13 @@ contains
real(8), intent(inout) :: buffer ! read data to here
logical, intent(in) :: collect ! collective I/O
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
if (collect) then
call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
else
call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, &
MPI_STATUS_IGNORE, mpiio_err)
end if
end subroutine mpi_read_double

51
src/output.F90
View file

@ -1372,7 +1372,7 @@ contains
if (cmfd_run) then
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
if (cmfd_display /= '') &
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========"
end if
write(UNIT=ou, FMT=*)
@ -1432,20 +1432,20 @@ contains
! write out cmfd keff if it is active and other display info
if (cmfd_on) then
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % k_cmfd(current_batch)
cmfd % k_cmfd(current_batch)
select case(trim(cmfd_display))
case('entropy')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % entropy(current_batch)
cmfd % entropy(current_batch)
case('balance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % balance(current_batch)
cmfd % balance(current_batch)
case('source')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % src_cmp(current_batch)
cmfd % src_cmp(current_batch)
case('dominance')
write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') &
cmfd % dom(current_batch)
cmfd % dom(current_batch)
end select
end if
@ -1567,7 +1567,7 @@ contains
speed_inactive = real(n_particles * (n_inactive - restart_batch) * &
gen_per_batch) / time_inactive % elapsed
speed_active = real(n_particles * n_active * gen_per_batch) / &
time_active % elapsed
time_active % elapsed
else
speed_inactive = ZERO
speed_active = real(n_particles * (n_batches - restart_batch) * &
@ -1884,21 +1884,22 @@ contains
select case(t % score_bins(k))
case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N)
score_name = 'P' // trim(to_str(t % moment_order(k))) // " " // &
score_names(abs(t % score_bins(k)))
score_names(abs(t % score_bins(k)))
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN)
score_index = score_index - 1
do n_order = 0, t % moment_order(k)
score_index = score_index + 1
score_name = 'P' // trim(to_str(n_order)) // " " //&
score_names(abs(t % score_bins(k)))
score_names(abs(t % score_bins(k)))
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) &
% sum_sq))
end do
k = k + t % moment_order(k)
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
@ -1908,11 +1909,13 @@ contains
do nm_order = -n_order, n_order
score_index = score_index + 1
score_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order)) // " " // score_names(abs(t % score_bins(k)))
trim(to_str(nm_order)) // " " &
// score_names(abs(t % score_bins(k)))
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index)&
% sum_sq))
end do
end do
k = k + (t % moment_order(k) + 1)**2 - 1
@ -1923,9 +1926,9 @@ contains
score_name = score_names(abs(t % score_bins(k)))
end if
write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') &
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
repeat(" ", indent), score_name, &
to_str(t % results(score_index,filter_index) % sum), &
trim(to_str(t % results(score_index,filter_index) % sum_sq))
end select
end do
indent = indent - 2
@ -2334,8 +2337,8 @@ contains
! Loop over lattice coordinates
do k = 1, n_x
do l = 1, n_y
do m = 1, n_z
do l = 1, n_y
do m = 1, n_z
if (final >= lat % offset(map, k, l, m) + offset) then
if (k == n_x .and. l == n_y .and. m == n_z) then

28
src/output_interface.F90
View file

@ -31,7 +31,7 @@ module output_interface
#endif
#endif
logical :: serial ! Serial I/O when using MPI/PHDF5
contains
contains
generic, public :: write_data => write_double, &
write_double_1Darray, &
write_double_2Darray, &
@ -111,7 +111,7 @@ contains
self % serial = serial
else
self % serial = .true.
end if
end if
#ifdef HDF5
# ifdef MPI
@ -153,7 +153,7 @@ contains
self % serial = serial
else
self % serial = .true.
end if
end if
#ifdef HDF5
# ifdef MPI
@ -201,18 +201,18 @@ contains
#ifdef HDF5
# ifdef MPI
call hdf5_file_close(self % hdf5_fh)
call hdf5_file_close(self % hdf5_fh)
# else
call hdf5_file_close(self % hdf5_fh)
call hdf5_file_close(self % hdf5_fh)
# endif
#elif MPI
if (self % serial) then
close(UNIT=self % unit_fh)
else
call mpi_close_file(self % mpi_fh)
end if
if (self % serial) then
close(UNIT=self % unit_fh)
else
call mpi_close_file(self % mpi_fh)
end if
#else
close(UNIT=self % unit_fh)
close(UNIT=self % unit_fh)
#endif
end subroutine file_close
@ -475,7 +475,7 @@ contains
call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length)
else
call hdf5_read_double_1Darray_parallel(self % hdf5_grp, name_, buffer, &
length, collect_)
length, collect_)
end if
# else
call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length)
@ -663,8 +663,8 @@ contains
if (self % serial) then
call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length)
else
call hdf5_write_double_3Darray_parallel(self % hdf5_grp, name_, buffer, length, &
collect_)
call hdf5_write_double_3Darray_parallel(self % hdf5_grp, name_, buffer, &
length, collect_)
end if
# else
call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length)

34
src/physics.F90
View file

@ -428,7 +428,7 @@ contains
! Sample velocity of target nucleus
if (.not. micro_xs(i_nuclide) % use_ptable) then
call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, &
& micro_xs(i_nuclide) % elastic)
& micro_xs(i_nuclide) % elastic)
else
v_t = ZERO
end if
@ -582,7 +582,7 @@ contains
! accompanying PDF and CDF is utilized)
if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. &
(sab % secondary_mode == SAB_SECONDARY_SKEWED)) then
(sab % secondary_mode == SAB_SECONDARY_SKEWED)) then
if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
! All bins equally likely
@ -843,16 +843,16 @@ contains
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
& maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
& maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up)
reject = .true.
@ -898,26 +898,26 @@ contains
! cdf value at lower bound attainable energy
if (i_E_low > 1) then
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
& + m * (E_low - nuc % energy_0K(i_E_low))
& + m * (E_low - nuc % energy_0K(i_E_low))
else
m = nuc % xs_cdf(i_E_low) &
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
& / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
end if
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
& / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
& + m * (E_up - nuc % energy_0K(i_E_up))
& + m * (E_up - nuc % energy_0K(i_E_up))
! values used to sample the Maxwellian
E_mode = kT
p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_mode / kT)
& * exp(-E_mode / kT)
E_t_max = 16.0_8 * E_mode
reject = .true.
@ -927,7 +927,7 @@ contains
! perform Maxwellian rejection sampling
E_t = E_t_max * prn()**2
p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) &
& * exp(-E_t / kT)
& * exp(-E_t / kT)
R_speed = p_t / p_mode
if (prn() < R_speed) then
@ -935,12 +935,12 @@ contains
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
& i_E_up - i_E_low + 2, cdf_rel)
& i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
& - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
& / (nuc % energy_0K(i_E_low + i_E_rel) &
& - nuc % energy_0K(i_E_low + i_E_rel - 1))
& - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
& / (nuc % energy_0K(i_E_low + i_E_rel) &
& - nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
! perform rejection sampling on cosine between
@ -1026,7 +1026,7 @@ contains
! Determine rejection probability
accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) &
/(beta_vn + beta_vt)
/(beta_vn + beta_vt)
! Perform rejection sampling on vt and mu
if (prn() < accept_prob) exit

2
src/plot.F90
View file

@ -250,7 +250,7 @@ contains
do j = ijk_ll(inner), ijk_ur(inner)
! check if we're in the mesh for this ijk
if (i > 0 .and. i <= m % dimension(outer) .and. &
j > 0 .and. j <= m % dimension(inner)) then
j > 0 .and. j <= m % dimension(inner)) then
! get xyz's of lower left and upper right of this mesh cell
xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1)

4
src/plot_header.F90
View file

@ -24,7 +24,7 @@ module plot_header
integer :: color_by ! quantity to color regions by
real(8) :: origin(3) ! xyz center of plot location
real(8) :: width(3) ! xyz widths of plot
integer :: basis ! direction of plot slice
integer :: basis ! direction of plot slice
integer :: pixels(3) ! pixel width/height of plot slice
integer :: meshlines_width ! pixel width of meshlines
integer :: level ! universe depth to plot the cells of
@ -37,7 +37,7 @@ module plot_header
! Plot type
integer, parameter :: PLOT_TYPE_SLICE = 1
integer, parameter :: PLOT_TYPE_VOXEL = 2
! Plot level
integer, parameter :: PLOT_LEVEL_LOWEST = -1

12
src/ppmlib.F90
View file

@ -1,7 +1,7 @@
module ppmlib
implicit none
!===============================================================================
! Image holds RGB information for output PPM image
!===============================================================================
@ -12,7 +12,7 @@ module ppmlib
end type Image
contains
!===============================================================================
! INIT_IMAGE initializes the Image derived type
!===============================================================================
@ -55,7 +55,7 @@ contains
!===============================================================================
! DEALLOCATE_IMAGE
!===============================================================================
subroutine deallocate_image(img)
type(Image) :: img
@ -70,7 +70,7 @@ contains
! INSIDE_IMAGE determines whether a point (x,y) is inside the image
!===============================================================================
function inside_image(img, x, y) result(inside)
type(Image), intent(in) :: img
@ -83,7 +83,7 @@ contains
(x >= 0) .and. (y >= 0)) inside = .true.
end function inside_image
!===============================================================================
! VALID_IMAGE checks whether the image has a width and height and if its color
! arrays are allocated
@ -123,5 +123,5 @@ contains
end if
end subroutine set_pixel
end module ppmlib

20
src/progress_header.F90
View file

@ -22,12 +22,12 @@ module progress_header
private
character(len=72) :: bar="???% | " // &
" |"
contains
procedure :: set_value => bar_set_value
contains
procedure :: set_value => bar_set_value
end type ProgressBar
contains
!===============================================================================
! IS_TERMINAL checks if output is currently being output to a terminal. Relies
! on a POSIX implementation of isatty, and defaults to false if that is not
@ -46,7 +46,7 @@ contains
#endif
end function is_terminal
!===============================================================================
! BAR_SET_VALUE prints the progress bar without advancing. The value is
! specified as percent completion, from 0 to 100. If the value is ever set to
@ -57,7 +57,7 @@ contains
class(ProgressBar), intent(inout) :: self
real(8), intent(in) :: val
integer :: i
if (.not. is_terminal()) return
@ -84,19 +84,19 @@ contains
write(OUTPUT_UNIT, '(A1,A1,A72)', ADVANCE='no') '+', char(13), self % bar
flush(OUTPUT_UNIT)
if (val >= 100.) then
! make new line
write(OUTPUT_UNIT, "(A)") ""
flush(OUTPUT_UNIT)
! reset the bar in case we want to use this instance again
self % bar = "???% | " // &
" |"
end if
end subroutine bar_set_value
end module progress_header

2
src/random_lcg.F90
View file

@ -93,7 +93,7 @@ contains
end do
end subroutine set_particle_seed
!===============================================================================
! PRN_SKIP advances the random number seed 'n' times from the current seed
!===============================================================================

6
src/search.F90
View file

@ -62,7 +62,7 @@ contains
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
&search.")
end if
end do
@ -114,7 +114,7 @@ contains
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
&search.")
end if
end do
@ -166,7 +166,7 @@ contains
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
&search.")
end if
end do

28
src/state_point.F90
View file

@ -52,7 +52,7 @@ contains
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
! Append appropriate extension
#ifdef HDF5
@ -256,7 +256,7 @@ contains
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j))
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT) then
tally % filters(j) % type == FILTER_ENERGYOUT) then
call sp % write_data(tally % filters(j) % real_bins, "bins", &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/filter " // to_str(j), &
@ -309,8 +309,8 @@ contains
do n_order = 0, tally % moment_order(k)
moment_name = 'P' // trim(to_str(n_order))
call sp % write_data(moment_name, "order" // trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
group="tallies/tally " // trim(to_str(tally % id)) // &
"/moments")
k = k + 1
end do
case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, &
@ -318,7 +318,7 @@ contains
do n_order = 0, tally % moment_order(k)
do nm_order = -n_order, n_order
moment_name = 'Y' // trim(to_str(n_order)) // ',' // &
trim(to_str(nm_order))
trim(to_str(nm_order))
call sp % write_data(moment_name, "order" // &
trim(to_str(k)), &
group="tallies/tally " // trim(to_str(tally % id)) // &
@ -412,7 +412,7 @@ contains
! Set filename
filename = trim(path_output) // 'source.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
#ifdef HDF5
filename = trim(filename) // '.h5'
@ -434,7 +434,7 @@ contains
! Set filename for state point
filename = trim(path_output) // 'statepoint.' // &
& zero_padded(current_batch, count_digits(n_max_batches))
& zero_padded(current_batch, count_digits(n_max_batches))
#ifdef HDF5
filename = trim(filename) // '.h5'
#else
@ -607,12 +607,12 @@ contains
tally % results(:,:) % sum_sq = tally_temp(2,:,:)
end if
! Put in temporary tally result
allocate(tallyresult_temp(m,n))
tallyresult_temp(:,:) % sum = tally_temp(1,:,:)
tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:)
! Put in temporary tally result
allocate(tallyresult_temp(m,n))
tallyresult_temp(:,:) % sum = tally_temp(1,:,:)
tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:)
! Write reduced tally results to file
! Write reduced tally results to file
call sp % write_tally_result(tally % results, "results", &
group="tallies/tally " // trim(to_str(tally % id)), n1=m, n2=n)
@ -687,7 +687,7 @@ contains
call sp % read_data(int_array(2), "version_minor")
call sp % read_data(int_array(3), "version_release")
if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR &
.or. int_array(3) /= VERSION_RELEASE) then
.or. int_array(3) /= VERSION_RELEASE) then
if (master) call warning("State point file was created with a different &
&version of OpenMC.")
end if
@ -843,7 +843,7 @@ contains
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j))
if (tally % filters(j) % type == FILTER_ENERGYIN .or. &
tally % filters(j) % type == FILTER_ENERGYOUT) then
tally % filters(j) % type == FILTER_ENERGYOUT) then
call sp % read_data(tally % filters(j) % real_bins, "bins", &
group="tallies/tally " // trim(to_str(curr_key)) // &
"/filter " // to_str(j), &

4
src/string.F90
View file

@ -7,7 +7,7 @@ module string
implicit none
interface to_str
module procedure int4_to_str, int8_to_str, real_to_str
module procedure int4_to_str, int8_to_str, real_to_str
end interface
contains
@ -50,7 +50,7 @@ contains
n = n + 1
if (i_end - i_start + 1 > len(words(n))) then
if (master) call warning("The word '" // string(i_start:i_end) &
&// "' is longer than the space allocated for it.")
&// "' is longer than the space allocated for it.")
end if
words(n) = string(i_start:i_end)
! reset indices

6
src/tally_header.F90
View file

@ -124,7 +124,7 @@ module tally_header
! Number of realizations of tally random variables
integer :: n_realizations = 0
! Tally precision triggers
integer :: n_triggers = 0 ! # of triggers
type(TriggerObject), allocatable :: triggers(:) ! Array of triggers
@ -197,10 +197,10 @@ module tally_header
this % reset = .false.
this % n_realizations = 0
if (allocated(this % triggers)) &
deallocate (this % triggers)
this % n_triggers = 0
end subroutine tallyobject_clear

12
src/timer_header.F90
View file

@ -15,11 +15,11 @@ module timer_header
logical :: running = .false. ! is timer running?
integer :: start_counts = 0 ! counts when started
real(8), public :: elapsed = ZERO ! total time elapsed in seconds
contains
procedure :: start => timer_start
procedure :: get_value => timer_get_value
procedure :: stop => timer_stop
procedure :: reset => timer_reset
contains
procedure :: start => timer_start
procedure :: get_value => timer_get_value
procedure :: stop => timer_stop
procedure :: reset => timer_reset
end type Timer
contains
@ -83,7 +83,7 @@ contains
!===============================================================================
subroutine timer_reset(self)
class(Timer), intent(inout) :: self
self % running = .false.

6
src/trigger.F90
View file

@ -43,7 +43,7 @@ contains
! When trigger threshold is reached, write information
if (satisfy_triggers) then
call write_message("Triggers satisfied for batch " // &
trim(to_str(current_batch)))
trim(to_str(current_batch)))
! When trigger is not reached write convergence info for user
elseif (name == "eigenvalue") then
@ -83,7 +83,7 @@ contains
! and finds the maximum uncertainty/threshold ratio for all triggers
!===============================================================================
subroutine check_tally_triggers(max_ratio, tally_id, name)
subroutine check_tally_triggers(max_ratio, tally_id, name)
! Variables to reflect distance to trigger convergence criteria
real(8), intent(inout) :: max_ratio ! max uncertainty/thresh ratio
@ -299,7 +299,7 @@ contains
! precision trigger(s).
!===============================================================================
subroutine compute_tally_current(t, trigger)
subroutine compute_tally_current(t, trigger)
integer :: i ! mesh index for x
integer :: j ! mesh index for y

6
src/trigger_header.F90
View file

@ -14,16 +14,16 @@ module trigger_header
character(len=52) :: score_name ! the name of the score
integer :: score_index ! the index of the score
real(8) :: variance=0.0 ! temp variance container
real(8) :: std_dev =0.0 ! temp std. dev. container
real(8) :: std_dev =0.0 ! temp std. dev. container
real(8) :: rel_err =0.0 ! temp rel. err. container
end type TriggerObject
!===============================================================================
! KTRIGGER describes a user-specified precision trigger for k-effective
!===============================================================================
type KTrigger
integer :: trigger_type = 0
real(8) :: threshold = 0
real(8) :: threshold = 0
end type KTrigger
end module trigger_header

12
src/vector_header.F90
View file

@ -9,12 +9,12 @@ module vector_header
integer :: n ! number of rows/cols in matrix
real(8), allocatable :: data(:) ! where vector data is stored
real(8), pointer :: val(:) ! pointer to vector data
contains
procedure :: create => vector_create
procedure :: destroy => vector_destroy
procedure :: add_value => vector_add_value
procedure :: copy => vector_copy
! TODO: procedure :: write => vector_write
contains
procedure :: create => vector_create
procedure :: destroy => vector_destroy
procedure :: add_value => vector_add_value
procedure :: copy => vector_copy
! TODO: procedure :: write => vector_write
end type Vector
contains

50
src/xml_interface.F90
View file

@ -92,7 +92,7 @@ contains
! Check if node exists
if (associated(temp_ptr)) return
! Check for a sub-element
! Check for a sub-element
elem_list => getChildrenByTagName(ptr, trim(node_name))
! Get the length of the list
@ -122,7 +122,7 @@ contains
! Set found to false
found_ = .false.
! Check for a sub-element
! Check for a sub-element
elem_list => getChildrenByTagName(in_ptr, trim(node_name))
! Get the length of the list
@ -147,7 +147,7 @@ contains
type(Node), pointer, intent(in) :: in_ptr
type(NodeList), pointer, intent(out) :: out_ptr
! Check for a sub-element
! Check for a sub-element
out_ptr => getChildrenByTagName(in_ptr, trim(node_name))
end subroutine get_node_list
@ -176,7 +176,7 @@ contains
type(NodeList), pointer, intent(in) :: in_ptr
type(Node), pointer, intent(out) :: out_ptr
! Check for a sub-element
! Check for a sub-element
out_ptr => item(in_ptr, idx - 1)
end subroutine get_list_item
@ -199,7 +199,7 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
@ -210,7 +210,7 @@ contains
else
call extractDataContent(temp_ptr, result_int)
end if
end subroutine get_node_value_integer
!===============================================================================
@ -231,18 +231,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_long)
else
call extractDataContent(temp_ptr, result_long)
end if
end subroutine get_node_value_long
!===============================================================================
@ -263,18 +263,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_double)
else
call extractDataContent(temp_ptr, result_double)
end if
end subroutine get_node_value_double
!===============================================================================
@ -295,18 +295,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_int)
else
call extractDataContent(temp_ptr, result_int)
end if
end subroutine get_node_array_integer
!===============================================================================
@ -327,18 +327,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_double)
else
call extractDataContent(temp_ptr, result_double)
end if
end subroutine get_node_array_double
!===============================================================================
@ -359,18 +359,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " &
&// getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_string)
else
call extractDataContent(temp_ptr, result_string)
end if
end subroutine get_node_array_string
!===============================================================================
@ -391,18 +391,18 @@ contains
call get_node(ptr, node_name, temp_ptr, node_type, found)
! Leave if it was not found
if (.not. found) then
if (.not. found) then
call fatal_error("Node " // node_name // " not part of Node " // &
getNodeName(ptr) // ".")
end if
! Extract value
if (node_type == ATTR_NODE) then
call extractDataAttribute(ptr, node_name, result_str)
else
call extractDataContent(temp_ptr, result_str)
end if
end subroutine get_node_value_string
!===============================================================================
@ -513,7 +513,7 @@ contains
! Check if node exists
if (associated(out_ptr)) return
! Check for a sub-element
! Check for a sub-element
elem_list => getChildrenByTagName(in_ptr, trim(node_name))
! Get the length of the list

232
tests/check_source.py Executable file
View file

@ -0,0 +1,232 @@
#!/usr/bin/env python
from __future__ import print_function
import glob
from string import whitespace
from sys import exit
from textwrap import TextWrapper
################################################################################
# Fortran reading/parsing backend.
################################################################################
class LineOfCode(object):
"""Contains and provides basic info about a line of Fortran 90 code."""
def __init__(self, number, content, is_continuation=False,
string_terminator=None):
# Initialize member variables.
self.number = number
self.content = content
self.is_continuation = is_continuation
self.is_continued = False
assert (string_terminator is None or
string_terminator == "'" or
string_terminator == '"')
self.initial_string_terminator = string_terminator
self.final_string_terminator = None
self.is_comment_only = False
# Parse the string.
self.parse()
def __repr__(self):
rep = 'LineOfCode: line number = {0:d}\n'.format(self.number)
rep += '\tis_continuation = ' + str(self.is_continuation) + '\n'
rep += '\tis_continued = ' + str(self.is_continued) + '\n'
rep += ('\tinitial_string_terminator = '
+ str(self.initial_string_terminator) + '\n')
rep += ('\tfinal_string_terminator = '
+ str(self.final_string_terminator) + '\n')
rep += '\tis_comment_only = ' + str(self.is_comment_only) + '\n'
rep += '\tContent:\n' + self.content
return rep
def __str__(self):
return repr(self)
def parse(self):
# Initialize the string variables.
if self.initial_string_terminator == "'":
in_a_single_quote = True
in_a_double_quote = False
elif self.initial_string_terminator == '"':
in_a_single_quote = False
in_a_double_quote = True
else:
in_a_single_quote = False
in_a_double_quote = False
# Initialize other variables.
in_a_comment = False
ends_with_ampersand = False
has_real_content = False
# Parse through the line.
for char in self.content:
# Check for the start or end of strings.
if char == "'" and in_a_single_quote:
in_a_single_quote = False
elif char == "'" and not in_a_double_quote:
in_a_single_quote = True
elif char == '"' and in_a_double_quote:
in_a_double_quote = False
elif char == '"' and not in_a_single_quote:
in_a_double_quote = True
# Check for the start of a comment.
if (char == "!" and not in_a_single_quote
and not in_a_double_quote):
in_a_comment = True
break # Don't care about anything in a comment.
# Check for a continuation ampersand.
if char == "&":
ends_with_ampersand = True
elif all([char != w for w in whitespace]):
ends_with_ampersand = False
# Check to see if there is any real content (non-whitespace, not in
# a comment)
if (not in_a_comment) and all([char != w for w in whitespace]):
has_real_content = True
# Make sure nothing went terribly wrong.
assert not (in_a_single_quote and in_a_double_quote)
assert not (in_a_single_quote and in_a_comment)
assert not (in_a_double_quote and in_a_comment)
# Is this line continued onto the next?
if ends_with_ampersand:
self.is_continued = True
# Is a multiline string continued on the next line?
if in_a_single_quote:
assert self.is_continued
self.final_string_terminator = "'"
if in_a_double_quote:
assert self.is_continued
self.final_string_terminator = '"'
# Is this line a comment-only line?
if (not has_real_content) and in_a_comment:
self.is_comment_only = True
def get_indent(self):
"""Return the number of indentation spaces prefixing this line."""
return len(self.content) - len(self.content.lstrip(' '))
def contains_whitespace_only(self):
"""Return True/False if all characters in the line are whitespace."""
if len(self.content.strip('\n')) == 0: return False
is_ws = [ any([char == ws for ws in whitespace])
for char in self.content.strip('\n') ]
return all(is_ws)
def has_trailing_whitespace(self):
"""Return True/False if this line ends with a whitespace character."""
stripped = self.content.strip('\n')
if len(stripped) == 0: return False
return any([stripped[-1] == ws for ws in whitespace])
def contains_tab(self):
"""Return True/False if a tab character appears in the line."""
return '\t' in self.content
def read_lines_of_code(fname):
line_num = 0
cont = False
str_term = None
with open(fname) as fin:
for line in fin:
loc = LineOfCode(line_num, line, is_continuation=cont,
string_terminator=str_term)
cont = loc.is_continued
str_term = loc.final_string_terminator
line_num += 1
yield loc
################################################################################
# Error checking.
################################################################################
def print_error(fname, line_number, err_msg):
header = "Error in file {0}, line {1}:".format(fname, line_number + 1)
tw = TextWrapper(width=80, initial_indent=' '*4, subsequent_indent=' '*4)
body = '\n'.join(tw.wrap(err_msg))
print(header + '\n' + body + '\n')
def check_source(fname):
"""Make sure the given Fortran source file meets OpenMC standards.
If errors are found, messages will be printed to the screen describing the
error. The function will return True if no errors were found or False
otherwise.
"""
good_code = True
base_indent = 0
for loc in read_lines_of_code(fname):
# Check for extra whitespace errors.
if loc.contains_whitespace_only():
good_code = False
print_error(fname, loc.number, "Line contains whitespace "
"characters but no content. Please remove whitespace.")
elif loc.has_trailing_whitespace():
good_code = False
print_error(fname, loc.number, "Line has trailing whitespace after"
" the content. Please remove trailing whitespace.")
if loc.contains_tab():
good_code = False
print_error(fname, loc.number, "Line contains a tab character. "
"Please replace with single whitespace characters.")
# Check indentation.
current_indent = loc.get_indent()
if ((not loc.is_continuation) and (not loc.is_comment_only) and
(not loc.contains_whitespace_only()) and current_indent % 2 != 0):
good_code = False
print_error(fname, loc.number, "Line is indented an odd number of "
"spaces. All non-continuation lines should be indented an "
"even number of spaces.")
if loc.is_continuation and current_indent < base_indent + 5:
good_code = False
print_error(fname, loc.number, "Continuation lines must be "
"indented by at least 5 spaces, but this line is indented {0}"
" spaces.".format(current_indent - base_indent))
# Set base indentation for next lines.
if not loc.is_continuation:
base_indent = current_indent
return good_code
################################################################################
# Main loop.
################################################################################
if __name__ == '__main__':
# Get a list of the F90 source files.
f90_files = glob.glob('../src/*.F90')
# Make sure we found the source files.
assert len(f90_files) != 0, 'No .F90 source files found.'
# Make sure all the F90 source files meet our standards.
good_code = [check_source(fname) for fname in f90_files]
if not all(good_code):
print("ERROR: The Fortran source code does not meet OpenMC's standards")
exit(-1)
else:
print("SUCCESS! Looks like the Fortran source meets our standards")
exit(0)

1
tests/travis.sh
View file

@ -3,6 +3,7 @@
set -ev
# Run all debug tests
./check_source.py
if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then
./run_tests.py -C "^basic-debug$|^hdf5-debug$|^mpi-omp-debug$|^phdf5-omp-debug$" -j 2 -s
else