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

Merged with upsream

Browse source
This commit is contained in:
Adam Nelson 2015-11-21 12:41:27 -05:00
commit d61e929dac
33 changed files with 1021 additions and 1130 deletions
Download patch file Download diff file Expand all files Collapse all files

62
openmc/mgxs/library.py
View file

@ -53,6 +53,8 @@ class Library(object):
The types of cross sections in the library (e.g., ['total', 'scatter'])
domain_type : {'material', 'cell', 'distribcell', 'universe'}
Domain type for spatial homogenization
domains : Iterable of Material, Cell or Universe
The spatial domain(s) for which MGXS in the Library are computed
correction : 'P0' or None
Apply the P0 correction to scattering matrices if set to 'P0'
energy_groups : EnergyGroups
@ -80,6 +82,7 @@ class Library(object):
self._by_nuclide = None
self._mgxs_types = []
self._domain_type = None
self._domains = 'all'
self._correction = 'P0'
self._energy_groups = None
self._tally_trigger = None
@ -105,6 +108,7 @@ class Library(object):
clone._by_nuclide = self.by_nuclide
clone._mgxs_types = self.mgxs_types
clone._domain_type = self.domain_type
clone._domains = self.domains
clone._correction = self.correction
clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
@ -153,22 +157,24 @@ class Library(object):
def by_nuclide(self):
return self._by_nuclide
@property
def domains(self):
if self.domain_type is None:
raise ValueError('Unable to get all domains without a domain type')
if self.domain_type == 'material':
return self.openmc_geometry.get_all_materials()
elif self.domain_type == 'cell' or self.domain_type == 'distribcell':
return self.openmc_geometry.get_all_material_cells()
elif self.domain_type == 'universe':
return self.openmc_geometry.get_all_universes()
@property
def domain_type(self):
return self._domain_type
@property
def domains(self):
if self._domains == 'all':
if self.domain_type == 'material':
return self.openmc_geometry.get_all_materials()
elif self.domain_type in ['cell', 'distribcell']:
return self.openmc_geometry.get_all_material_cells()
elif self.domain_type == 'universe':
return self.openmc_geometry.get_all_universes()
else:
raise ValueError('Unable to get domains without a domain type')
else:
return self._domains
@property
def correction(self):
return self._correction
@ -223,6 +229,38 @@ class Library(object):
cv.check_value('domain type', domain_type, tuple(openmc.mgxs.DOMAIN_TYPES))
self._domain_type = domain_type
@domains.setter
def domains(self, domains):
# Use all materials, cells or universes in the geometry as domains
if domains == 'all':
self._domains = domains
# User specified a list of material, cell or universe domains
else:
if self.domain_type == 'material':
cv.check_iterable_type('domain', domains, openmc.Material)
all_domains = self.openmc_geometry.get_all_materials()
elif self.domain_type in ['cell', 'distribcell']:
cv.check_iterable_type('domain', domains, openmc.Cell)
all_domains = self.openmc_geometry.get_all_material_cells()
elif self.domain_type == 'universe':
cv.check_iterable_type('domain', domains, openmc.Universe)
all_domains = self.openmc_geometry.get_all_universes()
else:
msg = 'Unable to set domains with ' \
'domain type "{}"'.format(self.domain_type)
raise ValueError(msg)
# Check that each domain can be found in the geometry
for domain in domains:
if domain not in all_domains:
msg = 'Domain "{}" could not be found in the ' \
'geometry.'.format(domain)
raise ValueError(msg)
self._domains = domains
@correction.setter
def correction(self, correction):
cv.check_value('correction', correction, ('P0', None))

34
openmc/opencg_compatible.py
View file

@ -486,20 +486,22 @@ def get_opencg_cell(openmc_cell):
# works if the region is a single half-space or an intersection of
# half-spaces, i.e., no complex cells.
region = openmc_cell.region
if isinstance(region, Halfspace):
surface = region.surface
halfspace = -1 if region.side == '-' else 1
opencg_cell.add_surface(get_opencg_surface(surface), halfspace)
elif isinstance(region, Intersection):
for node in region.nodes:
if not isinstance(node, Halfspace):
raise NotImplementedError("Complex cells not yet supported "
"in OpenCG.")
surface = node.surface
halfspace = -1 if node.side == '-' else 1
if region is not None:
if isinstance(region, Halfspace):
surface = region.surface
halfspace = -1 if region.side == '-' else 1
opencg_cell.add_surface(get_opencg_surface(surface), halfspace)
else:
raise NotImplementedError("Complex cells not yet supported in OpenCG.")
elif isinstance(region, Intersection):
for node in region.nodes:
if not isinstance(node, Halfspace):
raise NotImplementedError("Complex cells not yet "
"supported in OpenCG.")
surface = node.surface
halfspace = -1 if node.side == '-' else 1
opencg_cell.add_surface(get_opencg_surface(surface), halfspace)
else:
raise NotImplementedError("Complex cells not yet supported "
"in OpenCG.")
# Add the OpenMC Cell to the global collection of all OpenMC Cells
OPENMC_CELLS[cell_id] = openmc_cell
@ -877,6 +879,7 @@ def get_opencg_lattice(openmc_lattice):
pitch = openmc_lattice.pitch
lower_left = openmc_lattice.lower_left
universes = openmc_lattice.universes
outer = openmc_lattice.outer
if len(pitch) == 2:
new_pitch = np.ones(3, dtype=np.float64)
@ -909,6 +912,8 @@ def get_opencg_lattice(openmc_lattice):
opencg_lattice.dimension = dimension
opencg_lattice.width = pitch
opencg_lattice.universes = universe_array
if outer is not None:
opencg_lattice.outside = get_opencg_universe(outer)
offset = np.array(lower_left, dtype=np.float64) - \
((np.array(pitch, dtype=np.float64) *
@ -955,6 +960,7 @@ def get_openmc_lattice(opencg_lattice):
width = opencg_lattice.width
offset = opencg_lattice.offset
universes = opencg_lattice.universes
outer = opencg_lattice.outside
# Initialize an empty array for the OpenMC nested Universes in this Lattice
universe_array = np.ndarray(tuple(np.array(dimension)),
@ -985,6 +991,8 @@ def get_openmc_lattice(opencg_lattice):
openmc_lattice.pitch = width
openmc_lattice.universes = universe_array
openmc_lattice.lower_left = lower_left
if outer is not None:
openmc_lattice.outer = get_openmc_universe(outer)
# Add the OpenMC Lattice to the global collection of all OpenMC Lattices
OPENMC_LATTICES[lattice_id] = openmc_lattice

196
openmc/plots.py
View file

@ -5,9 +5,9 @@ import sys
import numpy as np
import openmc
import openmc.checkvalue as cv
from openmc.clean_xml import *
from openmc.checkvalue import (check_type, check_value, check_length,
check_greater_than, check_less_than)
if sys.version_info[0] >= 3:
basestring = str
@ -143,70 +143,70 @@ class Plot(object):
self._id = AUTO_PLOT_ID
AUTO_PLOT_ID += 1
else:
check_type('plot ID', plot_id, Integral)
check_greater_than('plot ID', plot_id, 0, equality=True)
cv.check_type('plot ID', plot_id, Integral)
cv.check_greater_than('plot ID', plot_id, 0, equality=True)
self._id = plot_id
@name.setter
def name(self, name):
check_type('plot name', name, basestring)
cv.check_type('plot name', name, basestring)
self._name = name
@width.setter
def width(self, width):
check_type('plot width', width, Iterable, Real)
check_length('plot width', width, 2, 3)
cv.check_type('plot width', width, Iterable, Real)
cv.check_length('plot width', width, 2, 3)
self._width = width
@origin.setter
def origin(self, origin):
check_type('plot origin', origin, Iterable, Real)
check_length('plot origin', origin, 3)
cv.check_type('plot origin', origin, Iterable, Real)
cv.check_length('plot origin', origin, 3)
self._origin = origin
@pixels.setter
def pixels(self, pixels):
check_type('plot pixels', pixels, Iterable, Integral)
check_length('plot pixels', pixels, 2, 3)
cv.check_type('plot pixels', pixels, Iterable, Integral)
cv.check_length('plot pixels', pixels, 2, 3)
for dim in pixels:
check_greater_than('plot pixels', dim, 0)
cv.check_greater_than('plot pixels', dim, 0)
self._pixels = pixels
@filename.setter
def filename(self, filename):
check_type('filename', filename, basestring)
cv.check_type('filename', filename, basestring)
self._filename = filename
@color.setter
def color(self, color):
check_type('plot color', color, basestring)
check_value('plot color', color, ['cell', 'mat'])
cv.check_type('plot color', color, basestring)
cv.check_value('plot color', color, ['cell', 'mat'])
self._color = color
@type.setter
def type(self, plottype):
check_type('plot type', plottype, basestring)
check_value('plot type', plottype, ['slice', 'voxel'])
cv.check_type('plot type', plottype, basestring)
cv.check_value('plot type', plottype, ['slice', 'voxel'])
self._type = plottype
@basis.setter
def basis(self, basis):
check_type('plot basis', basis, basestring)
check_value('plot basis', basis, ['xy', 'xz', 'yz'])
cv.check_type('plot basis', basis, basestring)
cv.check_value('plot basis', basis, ['xy', 'xz', 'yz'])
self._basis = basis
@background.setter
def background(self, background):
check_type('plot background', background, Iterable, Integral)
check_length('plot background', background, 3)
cv.check_type('plot background', background, Iterable, Integral)
cv.check_length('plot background', background, 3)
for rgb in background:
check_greater_than('plot background',rgb, 0, True)
check_less_than('plot background', rgb, 256)
cv.check_greater_than('plot background',rgb, 0, True)
cv.check_less_than('plot background', rgb, 256)
self._background = background
@col_spec.setter
def col_spec(self, col_spec):
check_type('plot col_spec parameter', col_spec, dict, Integral)
cv.check_type('plot col_spec parameter', col_spec, dict, Integral)
for key in col_spec:
if key < 0:
@ -229,18 +229,18 @@ class Plot(object):
@mask_componenets.setter
def mask_components(self, mask_components):
check_type('plot mask_components', mask_components, Iterable, Integral)
cv.check_type('plot mask_components', mask_components, Iterable, Integral)
for component in mask_components:
check_greater_than('plot mask_components', component, 0, True)
cv.check_greater_than('plot mask_components', component, 0, True)
self._mask_components = mask_components
@mask_background.setter
def mask_background(self, mask_background):
check_type('plot mask background', mask_background, Iterable, Integral)
check_length('plot mask background', mask_background, 3)
cv.check_type('plot mask background', mask_background, Iterable, Integral)
cv.check_length('plot mask background', mask_background, 3)
for rgb in mask_background:
check_greater_than('plot mask background', rgb, 0, True)
check_less_than('plot mask background', rgb, 256)
cv.check_greater_than('plot mask background', rgb, 0, True)
cv.check_less_than('plot mask background', rgb, 256)
self._mask_background = mask_background
def __repr__(self):
@ -261,6 +261,97 @@ class Plot(object):
string += '{0: <16}{1}{2}\n'.format('\tCol Spec', '=\t', self._col_spec)
return string
def colorize(self, geometry, seed=1):
"""Generate a color scheme for each domain in the plot.
This routine may be used to generate random, reproducible color schemes.
The colors generated are based upon cell/material IDs in the geometry.
Parameters
----------
geometry : openmc.Geometry
The geometry for which the plot is defined
seed : Integral
The random number seed used to generate the color scheme
"""
cv.check_type('geometry', geometry, openmc.Geometry)
cv.check_type('seed', seed, Integral)
cv.check_greater_than('seed', seed, 1, equality=True)
# Get collections of the domains which will be plotted
if self.color is 'mat':
domains = geometry.get_all_materials()
else:
domains = geometry.get_all_cells()
# Set the seed for the random number generator
np.random.seed(seed)
# Generate random colors for each feature
self.col_spec = {}
for domain in domains:
r = np.random.randint(0, 256)
g = np.random.randint(0, 256)
b = np.random.randint(0, 256)
self.col_spec[domain] = (r, g, b)
def highlight_domains(self, geometry, domains, seed=1,
alpha=0.5, background='gray'):
"""Use alpha compositing to highlight one or more domains in the plot.
This routine generates a color scheme and applies alpha compositing
to make all domains except the highlighted ones appear partially
transparent.
Parameters
----------
geometry : openmc.Geometry
The geometry for which the plot is defined
domains : Iterable of Integral
A collection of the domain IDs to highlight in the plot
seed : Integral
The random number seed used to generate the color scheme
alpha : Real in [0,1]
The value to apply in alpha compisiting
background : 3-tuple of Integral or 'white' or 'black' or 'gray'
The background color to apply in alpha compisiting
"""
cv.check_iterable_type('domains', domains, Integral)
cv.check_type('alpha', alpha, Real)
cv.check_greater_than('alpha', alpha, 0., equality=True)
cv.check_less_than('alpha', alpha, 1., equality=True)
# Get a background (R,G,B) tuple to apply in alpha compositing
if isinstance(background, basestring):
if background == 'white':
background = (255, 255, 255)
elif background == 'black':
background = (0, 0, 0)
elif background == 'gray':
background = (160, 160, 160)
else:
msg = 'The background "{}" is not defined'.format(background)
raise ValueError(msg)
cv.check_iterable_type('background', background, Integral)
# Generate a color scheme
self.colorize(geometry, seed)
# Apply alpha compositing to the colors for all domains
# other than those the user wishes to highlight
for domain_id in self.col_spec:
if domain_id not in domains:
r, g, b = self.col_spec[domain_id]
r = int(((1-alpha) * background[0]) + (alpha * r))
g = int(((1-alpha) * background[1]) + (alpha * g))
b = int(((1-alpha) * background[2]) + (alpha * b))
self._col_spec[domain_id] = (r, g, b)
def get_plot_xml(self):
"""Return XML representation of the plot
@ -349,6 +440,51 @@ class PlotsFile(object):
self._plots.remove(plot)
def colorize(self, geometry, seed=1):
"""Generate a consistent color scheme for each domain in each plot.
This routine may be used to generate random, reproducible color schemes.
The colors generated are based upon cell/material IDs in the geometry.
The color schemes will be consistent for all plots in "plots.xml".
Parameters
----------
geometry : openmc.Geometry
The geometry for which the plots are defined
seed : Integral
The random number seed used to generate the color scheme
"""
for plot in self._plots:
plot.colorize(geometry, seed)
def highlight_domains(self, geometry, domains, seed=1,
alpha=0.5, background='gray'):
"""Use alpha compositing to highlight one or more domains in the plot.
This routine generates a color scheme and applies alpha compositing
to make all domains except the highlighted ones partially transparent.
Parameters
----------
geometry : openmc.Geometry
The geometry for which the plot is defined
domains : Iterable of Integral
A collection of the domain IDs to highlight in the plot
seed : Integral
The random number seed used to generate the color scheme
alpha : Real in [0,1]
The value to apply in alpha compisiting
background : 3-tuple of Integral or 'white' or 'black' or 'gray'
The background color to apply in alpha compisiting
"""
for plot in self._plots:
plot.highlight_domains(geometry, domains, seed, alpha, background)
def _create_plot_subelements(self):
for plot in self._plots:
xml_element = plot.get_plot_xml()

86
openmc/settings.py
View file

@ -20,6 +20,8 @@ class SettingsFile(object):
Attributes
----------
run_mode : {'eigenvalue' or 'fixed source'}
The type of calculation to perform (default is 'eigenvalue')
batches : int
Number of batches to simulate
generations_per_batch : int
@ -128,7 +130,9 @@ class SettingsFile(object):
"""
def __init__(self):
# Eigenvalue subelement
# Run mode subelement (default is 'eigenvalue')
self._run_mode = 'eigenvalue'
self._batches = None
self._generations_per_batch = None
self._inactive = None
@ -206,9 +210,13 @@ class SettingsFile(object):
self._dd_count_interactions = False
self._settings_file = ET.Element("settings")
self._eigenvalue_subelement = None
self._run_mode_subelement = None
self._source_element = None
@property
def run_mode(self):
return self._run_mode
@property
def batches(self):
return self._batches
@ -417,6 +425,14 @@ class SettingsFile(object):
def dd_count_interactions(self):
return self._dd_count_interactions
@run_mode.setter
def run_mode(self, run_mode):
if 'run_mode' not in ['eigenvalue', 'fixed source']:
msg = 'Unable to set run mode to "{0}". Only "eigenvalue" ' \
'and "fixed source" are supported."'.format(run_mode)
raise ValueError(msg)
self._run_mode = run_mode
@batches.setter
def batches(self, batches):
check_type('batches', batches, Integral)
@ -891,72 +907,52 @@ class SettingsFile(object):
self._dd_count_interactions = interactions
def _create_eigenvalue_subelement(self):
self._create_particles_subelement()
self._create_batches_subelement()
self._create_inactive_subelement()
self._create_generations_per_batch_subelement()
self._create_keff_trigger_subelement()
def _create_run_mode_subelement(self):
if self.run_mode == 'eigenvalue':
self._run_mode_subelement = \
ET.SubElement(self._settings_file, "eigenvalue")
self._create_particles_subelement()
self._create_batches_subelement()
self._create_inactive_subelement()
self._create_generations_per_batch_subelement()
self._create_keff_trigger_subelement()
else:
if self._run_mode_subelement is None:
self._run_mode_subelement = \
ET.SubElement(self._settings_file, "fixed_source")
self._create_particles_subelement()
self._create_batches_subelement()
def _create_batches_subelement(self):
if self._batches is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "batches")
element = ET.SubElement(self._run_mode_subelement, "batches")
element.text = str(self._batches)
def _create_generations_per_batch_subelement(self):
if self._generations_per_batch is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement,
element = ET.SubElement(self._run_mode_subelement,
"generations_per_batch")
element.text = str(self._generations_per_batch)
def _create_inactive_subelement(self):
if self._inactive is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "inactive")
element = ET.SubElement(self._run_mode_subelement, "inactive")
element.text = str(self._inactive)
def _create_particles_subelement(self):
if self._particles is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "particles")
element = ET.SubElement(self._run_mode_subelement, "particles")
element.text = str(self._particles)
def _create_keff_trigger_subelement(self):
if self._keff_trigger is not None:
if self._eigenvalue_subelement is None:
self._eigenvalue_subelement = ET.SubElement(self._settings_file,
"eigenvalue")
element = ET.SubElement(self._eigenvalue_subelement, "keff_trigger")
element = ET.SubElement(self._run_mode_subelement, "keff_trigger")
for key in self._keff_trigger:
subelement = ET.SubElement(element, key)
subelement.text = str(self._keff_trigger[key]).lower()
def _create_energy_mode_subelement(self):
if self._energy_mode is not None:
element = ET.SubElement(self._settings_file, "energy_mode")
element.text = str(self._energy_mode)
def _create_max_order_subelement(self):
if self._max_order is not None:
element = ET.SubElement(self._settings_file, "max_order")
element.text = str(self._max_order)
def _create_source_subelement(self):
self._create_source_space_subelement()
self._create_source_energy_subelement()
@ -1222,10 +1218,10 @@ class SettingsFile(object):
self._settings_file.clear()
self._source_subelement = None
self._trigger_subelement = None
self._eigenvalue_subelement = None
self._run_mode_subelement = None
self._source_element = None
self._create_eigenvalue_subelement()
self._create_run_mode_subelement()
self._create_energy_mode_subelement()
self._create_max_order_subelement()
self._create_source_subelement()

404
src/ace.F90
View file

@ -46,9 +46,9 @@ contains
integer :: temp_table ! temporary value for sorting
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
type(Material), pointer :: mat => null()
type(Nuclide_CE), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(Material), pointer :: mat
type(Nuclide_CE), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(SetChar) :: already_read
! allocate arrays for ACE table storage and cross section cache
@ -235,7 +235,6 @@ contains
!===============================================================================
subroutine read_ace_table(i_table, i_listing)
integer, intent(in) :: i_table ! index in nuclides/sab_tables
integer, intent(in) :: i_listing ! index in xs_listings
@ -259,9 +258,9 @@ contains
character(10) :: mat ! material identifier
character(70) :: comment ! comment for ACE table
character(MAX_FILE_LEN) :: filename ! path to ACE cross section library
type(Nuclide_CE), pointer :: nuc => null()
type(SAlphaBeta), pointer :: sab => null()
type(XsListing), pointer :: listing => null()
type(Nuclide_CE), pointer :: nuc
type(SAlphaBeta), pointer :: sab
type(XsListing), pointer :: listing
! determine path, record length, and location of table
listing => xs_listings(i_listing)
@ -407,8 +406,6 @@ contains
end select
deallocate(XSS)
if(associated(nuc)) nullify(nuc)
if(associated(sab)) nullify(sab)
end subroutine read_ace_table
@ -418,10 +415,8 @@ contains
!===============================================================================
subroutine read_esz(nuc, data_0K)
type(Nuclide_CE), pointer :: nuc
logical :: data_0K ! are we reading 0K data?
type(Nuclide_CE), intent(inout) :: nuc
logical, intent(in) :: data_0K ! are we reading 0K data?
integer :: NE ! number of energy points for total and elastic cross sections
integer :: i ! index in 0K elastic xs array for this nuclide
@ -508,8 +503,7 @@ contains
!===============================================================================
subroutine read_nu_data(nuc)
type(Nuclide_CE), pointer :: nuc
type(Nuclide_CE), intent(inout) :: nuc
integer :: i ! loop index
integer :: JXS2 ! location for fission nu data
@ -525,7 +519,7 @@ contains
integer :: LOCC ! location of energy distributions for given MT
integer :: lc ! locator
integer :: length ! length of data to allocate
type(DistEnergy), pointer :: edist => null()
type(DistEnergy), pointer :: edist
JXS2 = JXS(2)
JXS24 = JXS(24)
@ -708,8 +702,7 @@ contains
!===============================================================================
subroutine read_reactions(nuc)
type(Nuclide_CE), pointer :: nuc
type(Nuclide_CE), intent(inout) :: nuc
integer :: i ! loop indices
integer :: i_fission ! index in nuc % index_fission
@ -723,7 +716,6 @@ contains
integer :: IE ! reaction's starting index on energy grid
integer :: NE ! number of energies
integer :: NR ! number of interpolation regions
type(Reaction), pointer :: rxn => null()
type(ListInt) :: MTs
LMT = JXS(3)
@ -741,14 +733,15 @@ contains
! Store elastic scattering cross-section on reaction one -- note that the
! sigma array is not allocated or stored for elastic scattering since it is
! already stored in nuc % elastic
rxn => nuc % reactions(1)
rxn % MT = 2
rxn % Q_value = ZERO
rxn % multiplicity = 1
rxn % threshold = 1
rxn % scatter_in_cm = .true.
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
associate (rxn => nuc % reactions(1))
rxn%MT = 2
rxn%Q_value = ZERO
rxn%multiplicity = 1
rxn%threshold = 1
rxn%scatter_in_cm = .true.
rxn%has_angle_dist = .false.
rxn%has_energy_dist = .false.
end associate
! Add contribution of elastic scattering to total cross section
nuc % total = nuc % total + nuc % elastic
@ -761,123 +754,125 @@ contains
i_fission = 0
do i = 1, NMT
rxn => nuc % reactions(i+1)
associate (rxn => nuc % reactions(i+1))
! set defaults
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
! set defaults
rxn % has_angle_dist = .false.
rxn % has_energy_dist = .false.
! read MT number, Q-value, and neutrons produced
rxn % MT = int(XSS(LMT + i - 1))
rxn % Q_value = XSS(JXS4 + i - 1)
rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1)))
rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0)
! read MT number, Q-value, and neutrons produced
rxn % MT = int(XSS(LMT + i - 1))
rxn % Q_value = XSS(JXS4 + i - 1)
rxn % multiplicity = abs(nint(XSS(JXS5 + i - 1)))
rxn % scatter_in_cm = (nint(XSS(JXS5 + i - 1)) < 0)
! Read energy-dependent multiplicities
if (rxn % multiplicity > 100) then
! Set flag and allocate space for Tab1 to store yield
rxn % multiplicity_with_E = .true.
allocate(rxn % multiplicity_E)
! Read energy-dependent multiplicities
if (rxn % multiplicity > 100) then
! 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
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
! allocate space for ENDF interpolation parameters
if (NR > 0) then
allocate(rxn % multiplicity_E % nbt(NR))
allocate(rxn % multiplicity_E % int(NR))
! 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)
rxn % multiplicity_E % y = get_real(NE)
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
! read starting energy index
LOCA = int(XSS(LXS + i - 1))
IE = int(XSS(JXS7 + LOCA - 1))
rxn % threshold = IE
! 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)
rxn % multiplicity_E % y = get_real(NE)
end if
! read starting energy index
LOCA = int(XSS(LXS + i - 1))
IE = int(XSS(JXS7 + LOCA - 1))
rxn % threshold = IE
! read number of energies cross section values
NE = int(XSS(JXS7 + LOCA))
allocate(rxn % sigma(NE))
XSS_index = JXS7 + LOCA + 1
rxn % sigma = get_real(NE)
! read number of energies cross section values
NE = int(XSS(JXS7 + LOCA))
allocate(rxn % sigma(NE))
XSS_index = JXS7 + LOCA + 1
rxn % sigma = get_real(NE)
end associate
end do
! Create set of MT values
do i = 1, size(nuc % reactions)
call MTs % append(nuc % reactions(i) % MT)
call nuc%reaction_index%add_key(nuc%reactions(i)%MT, i)
end do
! Create total, absorption, and fission cross sections
do i = 2, size(nuc % reactions)
rxn => nuc % reactions(i)
IE = rxn % threshold
NE = size(rxn % sigma)
associate (rxn => nuc % reactions(i))
IE = rxn % threshold
NE = size(rxn % sigma)
! Skip total inelastic level scattering, gas production cross sections
! (MT=200+), etc.
if (rxn % MT == N_LEVEL) cycle
if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle
! Skip total inelastic level scattering, gas production cross sections
! (MT=200+), etc.
if (rxn % MT == N_LEVEL) cycle
if (rxn % MT > N_5N2P .and. rxn % MT < N_P0) cycle
! Skip level cross sections if total is available
if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle
if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle
if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle
if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle
if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle
if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle
! Skip level cross sections if total is available
if (rxn % MT >= N_P0 .and. rxn % MT <= N_PC .and. MTs % contains(N_P)) cycle
if (rxn % MT >= N_D0 .and. rxn % MT <= N_DC .and. MTs % contains(N_D)) cycle
if (rxn % MT >= N_T0 .and. rxn % MT <= N_TC .and. MTs % contains(N_T)) cycle
if (rxn % MT >= N_3HE0 .and. rxn % MT <= N_3HEC .and. MTs % contains(N_3HE)) cycle
if (rxn % MT >= N_A0 .and. rxn % MT <= N_AC .and. MTs % contains(N_A)) cycle
if (rxn % MT >= N_2N0 .and. rxn % MT <= N_2NC .and. MTs % contains(N_2N)) cycle
! Add contribution to total cross section
nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma
! Add contribution to total cross section
nuc % total(IE:IE+NE-1) = nuc % total(IE:IE+NE-1) + rxn % sigma
! Add contribution to absorption cross section
if (is_disappearance(rxn % MT)) then
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
end if
! Add contribution to absorption cross section
if (is_disappearance(rxn % MT)) then
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
end if
! Information about fission reactions
if (rxn % MT == N_FISSION) then
allocate(nuc % index_fission(1))
elseif (rxn % MT == N_F) then
allocate(nuc % index_fission(PARTIAL_FISSION_MAX))
nuc % has_partial_fission = .true.
end if
! Information about fission reactions
if (rxn % MT == N_FISSION) then
allocate(nuc % index_fission(1))
elseif (rxn % MT == N_F) then
allocate(nuc % index_fission(PARTIAL_FISSION_MAX))
nuc % has_partial_fission = .true.
end if
! Add contribution to fission cross section
if (is_fission(rxn % MT)) then
nuc % fissionable = .true.
nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma
! Add contribution to fission cross section
if (is_fission(rxn % MT)) then
nuc % fissionable = .true.
nuc % fission(IE:IE+NE-1) = nuc % fission(IE:IE+NE-1) + rxn % sigma
! Also need to add fission cross sections to absorption
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
! Also need to add fission cross sections to absorption
nuc % absorption(IE:IE+NE-1) = nuc % absorption(IE:IE+NE-1) + rxn % sigma
! If total fission reaction is present, there's no need to store the
! reaction cross-section since it was copied to nuc % fission
if (rxn % MT == N_FISSION) deallocate(rxn % sigma)
! If total fission reaction is present, there's no need to store the
! reaction cross-section since it was copied to nuc % fission
if (rxn % MT == N_FISSION) deallocate(rxn % sigma)
! Keep track of this reaction for easy searching later
i_fission = i_fission + 1
nuc % index_fission(i_fission) = i
nuc % n_fission = nuc % n_fission + 1
end if
! Keep track of this reaction for easy searching later
i_fission = i_fission + 1
nuc % index_fission(i_fission) = i
nuc % n_fission = nuc % n_fission + 1
end if
end associate
end do
! Clear MTs set
@ -891,8 +886,7 @@ contains
!===============================================================================
subroutine read_angular_dist(nuc)
type(Nuclide_CE), pointer :: nuc
type(Nuclide_CE), intent(inout) :: nuc
integer :: JXS8 ! location of angular distribution locators
integer :: JXS9 ! location of angular distributions
@ -903,7 +897,6 @@ contains
integer :: i ! index in reactions array
integer :: j ! index over incoming energies
integer :: length ! length of data array to allocate
type(Reaction), pointer :: rxn => null()
JXS8 = JXS(8)
JXS9 = JXS(9)
@ -911,71 +904,72 @@ contains
! loop over all reactions with secondary neutrons -- NXS(5) does not include
! elastic scattering
do i = 1, NXS(5) + 1
rxn => nuc%reactions(i)
associate (rxn => nuc%reactions(i))
! find location of angular distribution
LOCB = int(XSS(JXS8 + i - 1))
if (LOCB == -1) then
! Angular distribution data are specified through LAWi = 44 in the DLW
! block
cycle
elseif (LOCB == 0) then
! No angular distribution data are given for this reaction, isotropic
! scattering is asssumed (in CM if TY < 0 and in LAB if TY > 0)
cycle
end if
rxn % has_angle_dist = .true.
! allocate space for incoming energies and locations
NE = int(XSS(JXS9 + LOCB - 1))
rxn % adist % n_energy = NE
allocate(rxn % adist % energy(NE))
allocate(rxn % adist % type(NE))
allocate(rxn % adist % location(NE))
! read incoming energy grid and location of nucs
XSS_index = JXS9 + LOCB
rxn % adist % energy = get_real(NE)
rxn % adist % location = get_int(NE)
! determine dize of data block
length = 0
do j = 1, NE
LC = rxn % adist % location(j)
if (LC == 0) then
! isotropic
rxn % adist % type(j) = ANGLE_ISOTROPIC
elseif (LC > 0) then
! 32 equiprobable bins
rxn % adist % type(j) = ANGLE_32_EQUI
length = length + 33
elseif (LC < 0) then
! tabular distribution
rxn % adist % type(j) = ANGLE_TABULAR
NP = int(XSS(JXS9 + abs(LC)))
length = length + 2 + 3*NP
! find location of angular distribution
LOCB = int(XSS(JXS8 + i - 1))
if (LOCB == -1) then
! Angular distribution data are specified through LAWi = 44 in the DLW
! block
cycle
elseif (LOCB == 0) then
! No angular distribution data are given for this reaction, isotropic
! scattering is assumed (in CM if TY < 0 and in LAB if TY > 0)
cycle
end if
end do
rxn % has_angle_dist = .true.
! allocate angular distribution data and read
allocate(rxn % adist % data(length))
! allocate space for incoming energies and locations
NE = int(XSS(JXS9 + LOCB - 1))
rxn % adist % n_energy = NE
allocate(rxn % adist % energy(NE))
allocate(rxn % adist % type(NE))
allocate(rxn % adist % location(NE))
! read angular distribution -- currently this does not actually parse the
! angular distribution tables for each incoming energy, that must be done
! on-the-fly
XSS_index = JXS9 + LOCB + 2 * NE
rxn % adist % data = get_real(length)
! read incoming energy grid and location of nucs
XSS_index = JXS9 + LOCB
rxn % adist % energy = get_real(NE)
rxn % adist % location = get_int(NE)
! change location pointers since they are currently relative to JXS(9)
LC = LOCB + 2 * NE + 1
do j = 1, NE
! For consistency, leave location as 0 if type is isotropic.
! This is not necessary for current correctness, but can avoid
! future issues
if (rxn % adist % location(j) /= 0) then
rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC
end if
end do
! determine dize of data block
length = 0
do j = 1, NE
LC = rxn % adist % location(j)
if (LC == 0) then
! isotropic
rxn % adist % type(j) = ANGLE_ISOTROPIC
elseif (LC > 0) then
! 32 equiprobable bins
rxn % adist % type(j) = ANGLE_32_EQUI
length = length + 33
elseif (LC < 0) then
! tabular distribution
rxn % adist % type(j) = ANGLE_TABULAR
NP = int(XSS(JXS9 + abs(LC)))
length = length + 2 + 3*NP
end if
end do
! allocate angular distribution data and read
allocate(rxn % adist % data(length))
! read angular distribution -- currently this does not actually parse the
! angular distribution tables for each incoming energy, that must be done
! on-the-fly
XSS_index = JXS9 + LOCB + 2 * NE
rxn % adist % data = get_real(length)
! change location pointers since they are currently relative to JXS(9)
LC = LOCB + 2 * NE + 1
do j = 1, NE
! For consistency, leave location as 0 if type is isotropic.
! This is not necessary for current correctness, but can avoid
! future issues
if (rxn % adist % location(j) /= 0) then
rxn % adist % location(j) = abs(rxn % adist % location(j)) - LC
end if
end do
end associate
end do
end subroutine read_angular_dist
@ -986,29 +980,28 @@ contains
!===============================================================================
subroutine read_energy_dist(nuc)
type(Nuclide_CE), pointer :: nuc
type(Nuclide_CE), intent(inout) :: nuc
integer :: LED ! location of energy distribution locators
integer :: LOCC ! location of energy distributions for given MT
integer :: i ! loop index
type(Reaction), pointer :: rxn => null()
LED = JXS(10)
! Loop over all reactions
do i = 1, NXS(5)
rxn => nuc % reactions(i+1) ! skip over elastic scattering
rxn % has_energy_dist = .true.
associate (rxn => nuc % reactions(i+1)) ! skip over elastic scattering
rxn % has_energy_dist = .true.
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
! find location of energy distribution data
LOCC = int(XSS(LED + i - 1))
! allocate energy distribution
allocate(rxn % edist)
! allocate energy distribution
allocate(rxn % edist)
! read data for energy distribution
call get_energy_dist(rxn % edist, LOCC)
! read data for energy distribution
call get_energy_dist(rxn % edist, LOCC)
end associate
end do
end subroutine read_energy_dist
@ -1020,10 +1013,9 @@ contains
!===============================================================================
recursive subroutine get_energy_dist(edist, loc_law, delayed_n)
type(DistEnergy), pointer :: edist ! energy distribution
integer, intent(in) :: loc_law ! locator for data
logical, optional :: delayed_n ! is this for delayed neutrons?
type(DistEnergy), intent(inout) :: edist ! energy distribution
integer, intent(in) :: loc_law ! locator for data
logical, intent(in), optional :: delayed_n ! is this for delayed neutrons?
integer :: LDIS ! location of all energy distributions
integer :: LNW ! location of next energy distribution if multiple
@ -1103,7 +1095,6 @@ contains
!===============================================================================
function length_energy_dist(lc, law, LOCC, lid) result(length)
integer, intent(in) :: lc ! location in XSS array
integer, intent(in) :: law ! energy distribution law
integer, intent(in) :: LOCC ! location of energy distribution
@ -1147,7 +1138,7 @@ contains
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
allocate(L(NE))
L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
@ -1205,7 +1196,7 @@ contains
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
allocate(L(NE))
L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
@ -1235,7 +1226,7 @@ contains
NR = int(XSS(lc + 1))
NE = int(XSS(lc + 2 + 2*NR))
allocate(L(NE))
L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Continue with finding data length
length = length + 2 + 2*NR + 2*NE
@ -1286,7 +1277,7 @@ contains
! in a way inconsistent with the current form of the ACE Format Guide
! (MCNP5 Manual, Vol 3)
allocate(L(NE))
L = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
L(:) = int(XSS(lc + 3 + 2*NR + NE: lc + 3 + 2*NR + 2*NE - 1))
! Don't currently do anything with L
deallocate(L)
! Continue with finding data length
@ -1302,8 +1293,7 @@ contains
!===============================================================================
subroutine read_unr_res(nuc)
type(Nuclide_CE), pointer :: nuc
type(Nuclide_CE), intent(inout) :: nuc
integer :: JXS23 ! location of URR data
integer :: lc ! locator
@ -1391,8 +1381,7 @@ contains
!===============================================================================
subroutine generate_nu_fission(nuc)
type(Nuclide_CE), pointer :: nuc
type(Nuclide_CE), intent(inout) :: nuc
integer :: i ! index on nuclide energy grid
real(8) :: E ! energy
@ -1418,8 +1407,7 @@ contains
!===============================================================================
subroutine read_thermal_data(table)
type(SAlphaBeta), pointer :: table
type(SAlphaBeta), intent(inout) :: table
integer :: i ! index for incoming energies
integer :: j ! index for outgoing energies
@ -1601,7 +1589,7 @@ contains
do i = 1, n_nuclides_total
do j = 1, n_nuclides_total
if (nuclides(i) % zaid == nuclides(j) % zaid) then
call nuclides(i) % nuc_list % append(j)
call nuclides(i) % nuc_list % push_back(j)
end if
end do
end do

56
src/ace_header.F90
View file

@ -1,8 +1,9 @@
module ace_header
use constants, only: MAX_FILE_LEN, ZERO
use endf_header, only: Tab1
use list_header, only: ListInt
use constants, only: MAX_FILE_LEN, ZERO
use dict_header, only: DictIntInt
use endf_header, only: Tab1
use stl_vector, only: VectorInt
implicit none
@ -17,10 +18,6 @@ module ace_header
integer, allocatable :: type(:) ! type of distribution
integer, allocatable :: location(:) ! location of each table
real(8), allocatable :: data(:) ! angular distribution data
! Type-Bound procedures
contains
procedure :: clear => distangle_clear ! Deallocates DistAngle
end type DistAngle
!===============================================================================
@ -51,7 +48,7 @@ module ace_header
integer :: MT ! ENDF MT value
real(8) :: Q_value ! Reaction Q value
integer :: multiplicity ! Number of secondary particles released
type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield
type(Tab1), allocatable :: multiplicity_E ! Energy-dependent neutron yield
integer :: threshold ! Energy grid index of threshold
logical :: scatter_in_cm ! scattering system in center-of-mass?
logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity
@ -79,27 +76,10 @@ module ace_header
logical :: multiply_smooth ! multiply by smooth cross section?
real(8), allocatable :: energy(:) ! incident energies
real(8), allocatable :: prob(:,:,:) ! actual probabibility tables
! Type-Bound procedures
contains
procedure :: clear => urrdata_clear ! Deallocates UrrData
end type UrrData
contains
!===============================================================================
! DISTANGLE_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine distangle_clear(this)
class(DistAngle), intent(inout) :: this ! The DistAngle object to clear
if (allocated(this % energy)) &
deallocate(this % energy, this % type, this % location, this % data)
end subroutine distangle_clear
!===============================================================================
! DISTENERGY_CLEAR resets and deallocates data in DistEnergy.
!===============================================================================
@ -108,12 +88,6 @@ module ace_header
class(DistEnergy), intent(inout) :: this ! The DistEnergy object to clear
! Clear p_valid
call this % p_valid % clear()
if (allocated(this % data)) &
deallocate(this % data)
if (associated(this % next)) then
! recursively clear this item
call this % next % clear()
@ -130,31 +104,11 @@ module ace_header
class(Reaction), intent(inout) :: this ! The Reaction object to clear
if (allocated(this % sigma)) deallocate(this % sigma)
if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E)
if (associated(this % edist)) then
call this % edist % clear()
deallocate(this % edist)
end if
call this % adist % clear()
end subroutine reaction_clear
!===============================================================================
! URRDATA_CLEAR resets and deallocates data in Reaction.
!===============================================================================
subroutine urrdata_clear(this)
class(UrrData), intent(inout) :: this ! The UrrData object to clear
if (allocated(this % energy)) &
deallocate(this % energy, this % prob)
end subroutine urrdata_clear
end module ace_header

103
src/cross_section.F90
View file

@ -15,10 +15,6 @@ module cross_section
use search, only: binary_search
implicit none
save
integer :: union_grid_index
!$omp threadprivate(union_grid_index)
contains
@ -29,13 +25,14 @@ contains
subroutine calculate_xs(p)
type(Particle), intent(in) :: p
type(Particle), intent(inout) :: p
integer :: i ! loop index over nuclides
integer :: i_nuclide ! index into nuclides array
integer :: i_sab ! index into sab_tables array
integer :: j ! index in mat % i_sab_nuclides
integer :: u ! index into logarithmic mapping array
integer :: i_grid ! index into logarithmic mapping array or material
! union grid
real(8) :: atom_density ! atom density of a nuclide
logical :: check_sab ! should we check for S(a,b) table?
type(Material), pointer :: mat ! current material
@ -46,7 +43,6 @@ contains
material_xs % absorption = ZERO
material_xs % fission = ZERO
material_xs % nu_fission = ZERO
material_xs % kappa_fission = ZERO
! Exit subroutine if material is void
if (p % material == MATERIAL_VOID) return
@ -54,11 +50,10 @@ contains
mat => materials(p % material)
! Find energy index on energy grid
u = 0
if (grid_method == GRID_MAT_UNION) then
call find_energy_index(p % E, p % material)
i_grid = find_energy_index(mat, p % E)
else if (grid_method == GRID_LOGARITHM) then
u = int(log(p % E/energy_min_neutron)/log_spacing)
i_grid = int(log(p % E/energy_min_neutron)/log_spacing)
end if
! Determine if this material has S(a,b) tables
@ -101,9 +96,9 @@ contains
! Calculate microscopic cross section for this nuclide
if (p % E /= micro_xs(i_nuclide) % last_E) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid)
else if (i_sab /= micro_xs(i_nuclide) % last_index_sab) then
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, u)
call calculate_nuclide_xs(i_nuclide, i_sab, p % E, p % material, i, i_grid)
end if
! ========================================================================
@ -131,10 +126,6 @@ contains
! Add contributions to material macroscopic nu-fission cross section
material_xs % nu_fission = material_xs % nu_fission + &
atom_density * micro_xs(i_nuclide) % nu_fission
! Add contributions to material macroscopic energy release from fission
material_xs % kappa_fission = material_xs % kappa_fission + &
atom_density * micro_xs(i_nuclide) % kappa_fission
end do
end subroutine calculate_xs
@ -144,20 +135,21 @@ contains
! given index in the nuclides array at the energy of the given particle
!===============================================================================
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, u)
subroutine calculate_nuclide_xs(i_nuclide, i_sab, E, i_mat, i_nuc_mat, i_log_union)
integer, intent(in) :: i_nuclide ! index into nuclides array
integer, intent(in) :: i_sab ! index into sab_tables array
real(8), intent(in) :: E ! energy
integer, intent(in) :: i_mat ! index into materials array
integer, intent(in) :: i_nuc_mat ! index into nuclides array for a material
integer, intent(in) :: u ! index into logarithmic mapping array
integer, intent(in) :: i_log_union ! index into logarithmic mapping array or
! material union energy grid
integer :: i_grid ! index on nuclide energy grid
integer :: i_low ! lower logarithmic mapping index
integer :: i_high ! upper logarithmic mapping index
real(8), intent(in) :: E ! energy
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide_CE), pointer :: nuc
type(Material), pointer :: mat
real(8) :: f ! interp factor on nuclide energy grid
type(Nuclide_CE), pointer :: nuc
type(Material), pointer :: mat
! Set pointer to nuclide and material
nuc => nuclides(i_nuclide)
@ -167,7 +159,7 @@ contains
select case (grid_method)
case (GRID_MAT_UNION)
i_grid = mat % nuclide_grid_index(i_nuc_mat, union_grid_index)
i_grid = mat % nuclide_grid_index(i_nuc_mat, i_log_union)
case (GRID_LOGARITHM)
! Determine the energy grid index using a logarithmic mapping to reduce
@ -180,8 +172,8 @@ contains
else
! Determine bounding indices based on which equal log-spaced interval
! the energy is in
i_low = nuc % grid_index(u)
i_high = nuc % grid_index(u + 1) + 1
i_low = nuc % grid_index(i_log_union)
i_high = nuc % grid_index(i_log_union + 1) + 1
! Perform binary search over reduced range
i_grid = binary_search(nuc % energy(i_low:i_high), &
@ -221,7 +213,6 @@ contains
! Initialize nuclide cross-sections to zero
micro_xs(i_nuclide) % fission = ZERO
micro_xs(i_nuclide) % nu_fission = ZERO
micro_xs(i_nuclide) % kappa_fission = ZERO
! Calculate microscopic nuclide total cross section
micro_xs(i_nuclide) % total = (ONE - f) * nuc % total(i_grid) &
@ -243,13 +234,6 @@ contains
! Calculate microscopic nuclide nu-fission cross section
micro_xs(i_nuclide) % nu_fission = (ONE - f) * nuc % nu_fission( &
i_grid) + f * nuc % nu_fission(i_grid+1)
! Calculate microscopic nuclide kappa-fission cross section
! The ENDF standard (ENDF-102) states that MT 18 stores
! the fission energy as the Q_value (fission(1))
micro_xs(i_nuclide) % kappa_fission = &
nuc % reactions(nuc % index_fission(1)) % Q_value * &
micro_xs(i_nuclide) % fission
end if
! If there is S(a,b) data for this nuclide, we need to do a few
@ -382,9 +366,8 @@ contains
real(8) :: fission ! fission cross section
real(8) :: inelastic ! inelastic cross section
logical :: same_nuc ! do we know the xs for this nuclide at this energy?
type(UrrData), pointer :: urr
type(UrrData), pointer :: urr
type(Nuclide_CE), pointer :: nuc
type(Reaction), pointer :: rxn
micro_xs(i_nuclide) % use_ptable = .true.
@ -410,7 +393,7 @@ contains
! preserve correlation of temperature in probability tables
same_nuc = .false.
do i = 1, nuc % nuc_list % size()
if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % get_item(i)) % last_E) then
if (E /= ZERO .and. E == micro_xs(nuc % nuc_list % data(i)) % last_E) then
same_nuc = .true.
same_nuc_idx = i
exit
@ -418,7 +401,7 @@ contains
end do
if (same_nuc) then
r = micro_xs(nuc % nuc_list % get_item(same_nuc_idx)) % last_prn
r = micro_xs(nuc % nuc_list % data(same_nuc_idx)) % last_prn
else
r = prn()
micro_xs(i_nuclide) % last_prn = r
@ -480,18 +463,17 @@ contains
! Determine treatment of inelastic scattering
inelastic = ZERO
if (urr % inelastic_flag > 0) then
! Get pointer to inelastic scattering reaction
rxn => nuc % reactions(nuc % urr_inelastic)
! Get index on energy grid and interpolation factor
i_energy = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
! Determine inelastic scattering cross section
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
end if
associate (rxn => nuc % reactions(nuc % urr_inelastic))
if (i_energy >= rxn % threshold) then
inelastic = (ONE - f) * rxn % sigma(i_energy - rxn%threshold + 1) + &
f * rxn % sigma(i_energy - rxn%threshold + 2)
end if
end associate
end if
! Multiply by smooth cross-section if needed
@ -528,38 +510,35 @@ contains
! energy
!===============================================================================
subroutine find_energy_index(E, i_mat)
real(8), intent(in) :: E ! energy of particle
integer, intent(in) :: i_mat ! material index
type(Material), pointer :: mat ! pointer to current material
mat => materials(i_mat)
pure function find_energy_index(mat, E) result(i)
type(Material), intent(in) :: mat ! pointer to current material
real(8), intent(in) :: E ! energy of particle
integer :: i ! energy grid index
! if the energy is outside of energy grid range, set to first or last
! index. Otherwise, do a binary search through the union energy grid.
if (E <= mat % e_grid(1)) then
union_grid_index = 1
i = 1
elseif (E > mat % e_grid(mat % n_grid)) then
union_grid_index = mat % n_grid - 1
i = mat % n_grid - 1
else
union_grid_index = binary_search(mat % e_grid, mat % n_grid, E)
i = binary_search(mat % e_grid, mat % n_grid, E)
end if
end subroutine find_energy_index
end function find_energy_index
!===============================================================================
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
! for a given nuclide at the trial relative energy used in resonance scattering
!===============================================================================
function elastic_xs_0K(E, nuc) result(xs_out)
pure function elastic_xs_0K(E, nuc) result(xs_out)
real(8), intent(in) :: E ! trial energy
type(Nuclide_CE), intent(in) :: nuc ! target nuclide at temperature
real(8) :: xs_out ! 0K xs at trial energy
type(Nuclide_CE), pointer :: nuc ! target nuclide at temperature
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
real(8), intent(inout) :: E ! trial energy
real(8) :: xs_out ! 0K xs at trial energy
integer :: i_grid ! index on nuclide energy grid
real(8) :: f ! interp factor on nuclide energy grid
! Determine index on nuclide energy grid
if (E < nuc % energy_0K(1)) then

2
src/endf.F90
View file

@ -11,7 +11,7 @@ contains
! REACTION_NAME gives the name of the reaction for a given MT value
!===============================================================================
function reaction_name(MT) result(string)
pure function reaction_name(MT) result(string)
integer, intent(in) :: MT
character(20) :: string

22
src/endf_header.F90
View file

@ -13,28 +13,6 @@ 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

30
src/fission.F90
View file

@ -15,10 +15,9 @@ contains
! given nuclide and incoming neutron energy
!===============================================================================
function nu_total(nuc, E) result(nu)
type(Nuclide_CE), pointer :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
pure function nu_total(nuc, E) result(nu)
type(Nuclide_CE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
integer :: i ! loop index
@ -26,7 +25,7 @@ contains
real(8) :: c ! polynomial coefficient
if (nuc % nu_t_type == NU_NONE) then
call fatal_error("No neutron emission data for table: " // nuc % name)
nu = ERROR_REAL
elseif (nuc % nu_t_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_t_data(1))
@ -49,10 +48,9 @@ contains
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_prompt(nuc, E) result(nu)
type(Nuclide_CE), pointer :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
pure function nu_prompt(nuc, E) result(nu)
type(Nuclide_CE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
integer :: i ! loop index
@ -87,11 +85,10 @@ contains
! for a given nuclide and incoming neutron energy
!===============================================================================
function nu_delayed(nuc, E) result(nu)
type(Nuclide_CE) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of delayed neutrons emitted per fission
pure function nu_delayed(nuc, E) result(nu)
type(Nuclide_CE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of delayed neutrons emitted per fission
if (nuc % nu_d_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron
@ -111,9 +108,8 @@ contains
! a given nuclide and incoming neutron energy in a given delayed group.
!===============================================================================
function yield_delayed(nuc, E, g) result(yield)
type(Nuclide_CE), pointer :: nuc ! nuclide from which to find nu
pure function yield_delayed(nuc, E, g) result(yield)
type(Nuclide_CE), intent(in) :: nuc ! nuclide from which to find nu
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: yield ! delayed neutron precursor yield
integer, intent(in) :: g ! the delayed neutron precursor group

111
src/geometry.F90
View file

@ -9,6 +9,7 @@ module geometry
use particle_header, only: LocalCoord, Particle
use particle_restart_write, only: write_particle_restart
use surface_header
use stl_vector, only: VectorInt
use simple_string, only: to_str
use tally, only: score_surface_current
@ -188,8 +189,8 @@ contains
recursive subroutine find_cell(p, found, search_cells)
type(Particle), intent(inout) :: p
logical, intent(inout) :: found
integer, optional :: search_cells(:)
logical, intent(inout) :: found
integer, optional :: search_cells(:)
integer :: i ! index over cells
integer :: j ! coordinate level index
integer :: i_xyz(3) ! indices in lattice
@ -199,7 +200,6 @@ contains
type(Cell), pointer :: c ! pointer to cell
class(Lattice), pointer :: lat ! pointer to lattice
type(Universe), pointer :: univ ! universe to search in
real(8) :: new_xyz(3) ! perturbed location used to look for cell
do j = p % n_coord + 1, MAX_COORD
call p % coord(j) % reset()
@ -286,8 +286,7 @@ contains
lat => lattices(c % fill) % obj
! Determine lattice indices
new_xyz = p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw
i_xyz = lat % get_indices(new_xyz)
i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw)
! Store lower level coordinates
p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz)
@ -342,7 +341,7 @@ contains
subroutine cross_surface(p, last_cell)
type(Particle), intent(inout) :: p
integer, intent(in) :: last_cell ! last cell particle was in
integer, intent(in) :: last_cell ! last cell particle was in
real(8) :: u ! x-component of direction
real(8) :: v ! y-component of direction
@ -504,7 +503,7 @@ contains
subroutine cross_lattice(p, lattice_translation)
type(Particle), intent(inout) :: p
integer, intent(in) :: lattice_translation(3)
integer, intent(in) :: lattice_translation(3)
integer :: j
integer :: i_xyz(3) ! indices in lattice
logical :: found ! particle found in cell?
@ -577,10 +576,10 @@ contains
subroutine distance_to_boundary(p, dist, surface_crossed, lattice_translation, &
next_level)
type(Particle), intent(inout) :: p
real(8), intent(out) :: dist
integer, intent(out) :: surface_crossed
integer, intent(out) :: lattice_translation(3)
integer, intent(out) :: next_level
real(8), intent(out) :: dist
integer, intent(out) :: surface_crossed
integer, intent(out) :: lattice_translation(3)
integer, intent(out) :: next_level
integer :: i ! index for surface in cell
integer :: j
@ -873,81 +872,51 @@ contains
subroutine neighbor_lists()
integer :: i ! index in cells/surfaces array
integer :: j ! index of surface in cell
integer :: i_surface ! index in count arrays
integer, allocatable :: count_positive(:) ! # of cells on positive side
integer, allocatable :: count_negative(:) ! # of cells on negative side
logical :: positive ! positive side specified in surface list
type(Cell), pointer :: c
integer :: i ! index in cells/surfaces array
integer :: j ! index in region specification
integer :: k ! surface half-space spec
integer :: n ! size of vector
type(VectorInt), allocatable :: neighbor_pos(:)
type(VectorInt), allocatable :: neighbor_neg(:)
call write_message("Building neighboring cells lists for each surface...", &
&4)
4)
allocate(count_positive(n_surfaces))
allocate(count_negative(n_surfaces))
count_positive = 0
count_negative = 0
allocate(neighbor_pos(n_surfaces))
allocate(neighbor_neg(n_surfaces))
do i = 1, n_cells
c => cells(i)
do j = 1, size(cells(i)%region)
! Get token from region specification and skip any tokens that
! correspond to operators rather than regions
k = cells(i)%region(j)
if (abs(k) >= OP_UNION) cycle
! loop over each region specification
do j = 1, size(c%region)
i_surface = c % region(j)
positive = (i_surface > 0)
! Skip any tokens that correspond to operators rather than regions
i_surface = abs(i_surface)
if (i_surface >= OP_UNION) cycle
if (positive) then
count_positive(i_surface) = count_positive(i_surface) + 1
! Add this cell ID to neighbor list for k-th surface
if (k > 0) then
call neighbor_pos(abs(k))%push_back(i)
else
count_negative(i_surface) = count_negative(i_surface) + 1
call neighbor_neg(abs(k))%push_back(i)
end if
end do
end do
! allocate neighbor lists for each surface
do i = 1, n_surfaces
if (count_positive(i) > 0) then
allocate(surfaces(i)%obj%neighbor_pos(count_positive(i)))
! Copy positive neighbors to Surface instance
n = neighbor_pos(i)%size()
if (n > 0) then
allocate(surfaces(i)%obj%neighbor_pos(n))
surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:n)
end if
if (count_negative(i) > 0) then
allocate(surfaces(i)%obj%neighbor_neg(count_negative(i)))
! Copy negative neighbors to Surface instance
n = neighbor_neg(i)%size()
if (n > 0) then
allocate(surfaces(i)%obj%neighbor_neg(n))
surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:n)
end if
end do
count_positive = 0
count_negative = 0
! loop over all cells
do i = 1, n_cells
c => cells(i)
! loop through the region specification
do j = 1, size(c%region)
i_surface = c % region(j)
positive = (i_surface > 0)
! Skip any tokens that correspond to operators rather than regions
i_surface = abs(i_surface)
if (i_surface >= OP_UNION) cycle
if (positive) then
count_positive(i_surface) = count_positive(i_surface) + 1
surfaces(i_surface)%obj%neighbor_pos(count_positive(i_surface)) = i
else
count_negative(i_surface) = count_negative(i_surface) + 1
surfaces(i_surface)%obj%neighbor_neg(count_negative(i_surface)) = i
end if
end do
end do
deallocate(count_positive)
deallocate(count_negative)
end subroutine neighbor_lists
!===============================================================================
@ -957,7 +926,7 @@ contains
subroutine handle_lost_particle(p, message)
type(Particle), intent(inout) :: p
character(*) :: message
character(*) :: message
! Print warning and write lost particle file
call warning(message)

16
src/global.F90
View file

@ -473,7 +473,12 @@ contains
do i = 1, size(nuclides)
call nuclides(i) % clear()
end do
deallocate(nuclides)
! WARNING: The following statement should work but doesn't under gfortran
! 4.6 because of a bug. Technically, commenting this out leaves a memory
! leak.
! deallocate(nuclides)
end if
if (allocated(nuclides_0K)) then
@ -516,14 +521,7 @@ contains
! Deallocate tally-related arrays
if (allocated(global_tallies)) deallocate(global_tallies)
if (allocated(meshes)) deallocate(meshes)
if (allocated(tallies)) then
! First call the clear routines
do i = 1, size(tallies)
call tallies(i) % clear()
end do
! Now deallocate the tally array
deallocate(tallies)
end if
if (allocated(tallies)) deallocate(tallies)
if (allocated(matching_bins)) deallocate(matching_bins)
if (allocated(tally_maps)) deallocate(tally_maps)

2
src/input_xml.F90
View file

@ -1037,7 +1037,7 @@ contains
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
character(MAX_LINE_LEN) :: region_spec
character(1000) :: region_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
class(Lattice), pointer :: lat

32
src/interpolation.F90
View file

@ -21,7 +21,7 @@ contains
! tabulated x's and y's.
!===============================================================================
function interpolate_tab1_array(data, x, loc_start) result(y)
pure function interpolate_tab1_array(data, x, loc_start) result(y)
real(8), intent(in) :: data(:) ! array of data
real(8), intent(in) :: x ! x value to find y at
@ -106,18 +106,16 @@ contains
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = (1 - r)*y0 + r*y1
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = (1 - r)*y0 + r*y1
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = exp((1-r)*log(y0) + r*log(y1))
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
call fatal_error("Unsupported interpolation scheme: " // to_str(interp))
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_array
@ -129,7 +127,7 @@ contains
! tabulated x's and y's.
!===============================================================================
function interpolate_tab1_object(obj, x) result(y)
pure function interpolate_tab1_object(obj, x) result(y)
type(Tab1), intent(in) :: obj ! ENDF Tab1 interpolable function
real(8), intent(in) :: x ! x value to find y at
@ -191,18 +189,16 @@ contains
select case (interp)
case (LINEAR_LINEAR)
r = (x - x0)/(x1 - x0)
y = (1 - r)*y0 + r*y1
y = y0 + r*(y1 - y0)
case (LINEAR_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = (1 - r)*y0 + r*y1
r = log(x/x0)/log(x1/x0)
y = y0 + r*(y1 - y0)
case (LOG_LINEAR)
r = (x - x0)/(x1 - x0)
y = exp((1-r)*log(y0) + r*log(y1))
y = y0*exp(r*log(y1/y0))
case (LOG_LOG)
r = (log(x) - log(x0))/(log(x1) - log(x0))
y = exp((1-r)*log(y0) + r*log(y1))
case default
call fatal_error("Unsupported interpolation scheme: " // to_str(interp))
r = log(x/x0)/log(x1/x0)
y = y0*exp(r*log(y1/y0))
end select
end function interpolate_tab1_object

6
src/macroxs.F90
View file

@ -34,9 +34,6 @@ contains
xs % absorption = this % absorption(gin)
xs % fission = this % fission(gin)
xs % nu_fission = this % nu_fission(gin)
if (allocated(this % k_fission)) then
xs % kappa_fission = this % k_fission(gin)
end if
type is (MacroXS_Angle)
call find_angle(this % polar, this % azimuthal, uvw, iazi, ipol)
@ -45,9 +42,6 @@ contains
xs % absorption = this % absorption(gin, iazi, ipol)
xs % fission = this % fission(gin, iazi, ipol)
xs % nu_fission = this % nu_fission(gin, iazi, ipol)
if (allocated(this % k_fission)) then
xs % kappa_fission = this % k_fission(gin, iazi, ipol)
end if
end select
end subroutine calculate_mgxs

6
src/math.F90
View file

@ -11,7 +11,7 @@ contains
! distribution with a specified probability level
!===============================================================================
function normal_percentile(p) result(z)
elemental function normal_percentile(p) result(z)
real(8), intent(in) :: p ! probability level
real(8) :: z ! corresponding z-value
@ -70,7 +70,7 @@ contains
! specified probability level and number of degrees of freedom
!===============================================================================
function t_percentile(p, df) result(t)
elemental function t_percentile(p, df) result(t)
real(8), intent(in) :: p ! probability level
integer, intent(in) :: df ! degrees of freedom
@ -122,7 +122,7 @@ contains
! the return value will be 1.0.
!===============================================================================
pure function calc_pn(n,x) result(pnx)
elemental function calc_pn(n,x) result(pnx)
integer, intent(in) :: n ! Legendre order requested
real(8), intent(in) :: x ! Independent variable the Legendre is to be

38
src/mesh.F90
View file

@ -17,9 +17,8 @@ contains
! GET_MESH_BIN determines the tally bin for a particle in a structured mesh
!===============================================================================
subroutine get_mesh_bin(m, xyz, bin)
type(RegularMesh), pointer :: m ! mesh pointer
pure subroutine get_mesh_bin(m, xyz, bin)
type(RegularMesh), intent(in) :: m ! mesh pointer
real(8), intent(in) :: xyz(:) ! coordinates
integer, intent(out) :: bin ! tally bin
@ -70,9 +69,8 @@ contains
! GET_MESH_INDICES determines the indices of a particle in a structured mesh
!===============================================================================
subroutine get_mesh_indices(m, xyz, ijk, in_mesh)
type(RegularMesh), pointer :: m
pure subroutine get_mesh_indices(m, xyz, ijk, in_mesh)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz(:) ! coordinates to check
integer, intent(out) :: ijk(:) ! indices in mesh
logical, intent(out) :: in_mesh ! were given coords in mesh?
@ -95,11 +93,10 @@ contains
! use in a TallyObject results array
!===============================================================================
function mesh_indices_to_bin(m, ijk, surface_current) result(bin)
type(RegularMesh), pointer :: m
pure function mesh_indices_to_bin(m, ijk, surface_current) result(bin)
type(RegularMesh), intent(in) :: m
integer, intent(in) :: ijk(:)
logical, optional :: surface_current
logical, intent(in), optional :: surface_current
integer :: bin
integer :: n_y ! number of mesh cells in y direction
@ -129,9 +126,8 @@ contains
! (i,j) or (i,j,k) indices
!===============================================================================
subroutine bin_to_mesh_indices(m, bin, ijk)
type(RegularMesh), pointer :: m
pure subroutine bin_to_mesh_indices(m, bin, ijk)
type(RegularMesh), intent(in) :: m
integer, intent(in) :: bin
integer, intent(out) :: ijk(:)
@ -166,9 +162,9 @@ contains
type(Bank), intent(in) :: bank_array(:) ! fission or source bank
real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each
! cell and energy group
real(8), optional :: energies(:) ! energy grid to search
integer(8), optional :: size_bank ! # of bank sites (on each proc)
logical, optional :: sites_outside ! were there sites outside mesh?
real(8), intent(in), optional :: energies(:) ! energy grid to search
integer(8), intent(in), optional :: size_bank ! # of bank sites (on each proc)
logical, intent(inout), optional :: sites_outside ! were there sites outside mesh?
integer :: i ! loop index for local fission sites
integer :: n_sites ! size of bank array
@ -261,9 +257,8 @@ contains
! track will score to a mesh tally.
!===============================================================================
function mesh_intersects_2d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), pointer :: m
pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(2)
real(8), intent(in) :: xyz1(2)
logical :: intersects
@ -327,9 +322,8 @@ contains
end function mesh_intersects_2d
function mesh_intersects_3d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), pointer :: m
pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects)
type(RegularMesh), intent(in) :: m
real(8), intent(in) :: xyz0(3)
real(8), intent(in) :: xyz1(3)
logical :: intersects

2
src/mgxs_data.F90
View file

@ -180,7 +180,7 @@ contains
do i = 1, n_nuclides_total
do j = 1, n_nuclides_total
if (nuclides_MG(i) % obj % zaid == nuclides_MG(j) % obj % zaid) then
call nuclides_MG(i) % obj % nuc_list % append(j)
call nuclides_MG(i) % obj % nuc_list % push_back(j)
end if
end do
end do

113
src/nuclide_header.F90
View file

@ -26,7 +26,7 @@ module nuclide_header
real(8) :: kT ! temperature in MeV (k*T)
! Linked list of indices in nuclides array of instances of this same nuclide
type(ListInt) :: nuc_list
type(VectorInt) :: nuc_list
! Fission information
logical :: fissionable ! nuclide is fissionable?
@ -37,6 +37,12 @@ module nuclide_header
end type Nuclide_Base
abstract interface
subroutine nuclide_base_clear_(this)
import Nuclide_Base
class(Nuclide_Base), intent(inout) :: this
end subroutine nuclide_base_clear_
subroutine print_nuclide_(this, unit)
import Nuclide_Base
class(Nuclide_Base),intent(in) :: this
@ -97,7 +103,9 @@ module nuclide_header
! Reactions
integer :: n_reaction ! # of reactions
type(Reaction), pointer :: reactions(:) => null()
type(Reaction), allocatable :: reactions(:)
type(DictIntInt) :: reaction_index ! map MT values to index in reactions
! array; used at tally-time
! Type-Bound procedures
contains
@ -243,7 +251,6 @@ module nuclide_header
real(8) :: absorption ! microscopic absorption xs
real(8) :: fission ! microscopic fission xs
real(8) :: nu_fission ! microscopic production xs
real(8) :: kappa_fission ! microscopic energy-released from fission
! Information for S(a,b) use
integer :: index_sab ! index in sab_tables (zero means no table)
@ -266,7 +273,6 @@ module nuclide_header
real(8) :: absorption ! macroscopic absorption xs
real(8) :: fission ! macroscopic fission xs
real(8) :: nu_fission ! macroscopic production xs
real(8) :: kappa_fission ! macroscopic energy-released from fission
end type MaterialMacroXS
!===============================================================================
@ -296,44 +302,12 @@ module nuclide_header
! or Nuclide_Angle
!===============================================================================
subroutine nuclide_base_clear_(this)
class(Nuclide_Base), intent(inout) :: this
call this % nuc_list % clear()
end subroutine nuclide_base_clear_
subroutine nuclide_ce_clear(this)
class(Nuclide_CE), intent(inout) :: this ! The Nuclide object to clear
integer :: i ! Loop counter
if (allocated(this % energy)) &
deallocate(this % energy, this % total, this % elastic, &
& this % fission, this % nu_fission, this % absorption)
if (allocated(this % energy_0K)) &
deallocate(this % energy_0K)
if (allocated(this % elastic_0K)) &
deallocate(this % elastic_0K)
if (allocated(this % xs_cdf)) &
deallocate(this % xs_cdf)
if (allocated(this % heating)) &
deallocate(this % heating)
if (allocated(this % index_fission)) deallocate(this % index_fission)
if (allocated(this % nu_t_data)) deallocate(this % nu_t_data)
if (allocated(this % nu_p_data)) deallocate(this % nu_p_data)
if (allocated(this % nu_d_data)) deallocate(this % nu_d_data)
if (allocated(this % nu_d_precursor_data)) &
deallocate(this % nu_d_precursor_data)
if (associated(this % nu_d_edist)) then
do i = 1, size(this % nu_d_edist)
call this % nu_d_edist(i) % clear()
@ -342,18 +316,16 @@ module nuclide_header
end if
if (associated(this % urr_data)) then
call this % urr_data % clear()
deallocate(this % urr_data)
end if
if (associated(this % reactions)) then
if (allocated(this % reactions)) then
do i = 1, size(this % reactions)
call this % reactions(i) % clear()
end do
deallocate(this % reactions)
end if
call nuclide_base_clear_(this)
call this % reaction_index % clear()
end subroutine nuclide_ce_clear
@ -361,9 +333,6 @@ module nuclide_header
class(Nuclide_Iso), intent(inout) :: this ! The Nuclide object to clear
! Clear the base object
call nuclide_base_clear_(this)
! Cler the extended information
if (allocated(this % total)) then
deallocate(this % total, this % absorption, this % scatter)
@ -387,9 +356,6 @@ module nuclide_header
class(Nuclide_Angle), intent(inout) :: this ! The Nuclide object to clear
! Clear the base object
call nuclide_base_clear_(this)
! Cler the extended information
if (allocated(this % total)) then
deallocate(this % total, this % absorption, this % scatter)
@ -436,8 +402,7 @@ module nuclide_header
integer :: size_energy ! memory used for a energy distributions (bytes)
integer :: size_urr ! memory used for probability tables (bytes)
character(11) :: law ! secondary energy distribution law
type(Reaction), pointer :: rxn => null()
type(UrrData), pointer :: urr => null()
type(UrrData), pointer :: urr
! set default unit for writing information
if (present(unit)) then
@ -465,32 +430,33 @@ module nuclide_header
! Information on each reaction
write(unit_,*) ' Reaction Q-value COM Law IE size(angle) size(energy)'
do i = 1, this % n_reaction
rxn => this % reactions(i)
associate (rxn => this % reactions(i))
! Determine size of angle distribution
if (rxn % has_angle_dist) then
size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
else
size_angle = 0
end if
! Determine size of angle distribution
if (rxn % has_angle_dist) then
size_angle = rxn % adist % n_energy * 16 + size(rxn % adist % data) * 8
else
size_angle = 0
end if
! Determine size of energy distribution and law
if (rxn % has_energy_dist) then
size_energy = size(rxn % edist % data) * 8
law = to_str(rxn % edist % law)
else
size_energy = 0
law = 'None'
end if
! Determine size of energy distribution and law
if (rxn % has_energy_dist) then
size_energy = size(rxn % edist % data) * 8
law = to_str(rxn % edist % law)
else
size_energy = 0
law = 'None'
end if
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,A4,1X,I6,1X,I11,1X,I11)') &
reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, &
law(1:4), rxn % threshold, size_angle, size_energy
! Accumulate data size
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
size_angle_total = size_angle_total + size_angle
size_energy_total = size_energy_total + size_energy
! Accumulate data size
size_xs = size_xs + (this % n_grid - rxn%threshold + 1) * 8
size_angle_total = size_angle_total + size_angle
size_energy_total = size_energy_total + size_energy
end associate
end do
! Add memory required for summary reactions (total, absorption, fission,
@ -574,7 +540,6 @@ module nuclide_header
integer :: unit_ ! unit to write to
integer :: size_total, size_scattmat, size_mgxs
character(MAX_LINE_LEN) :: temp_str
! set default unit for writing information
if (present(unit)) then
@ -615,7 +580,6 @@ module nuclide_header
integer :: unit_ ! unit to write to
integer :: size_total, size_scattmat, size_mgxs
character(MAX_LINE_LEN) :: temp_str
! set default unit for writing information
if (present(unit)) then
@ -668,9 +632,6 @@ module nuclide_header
integer, optional, intent(in) :: i_pol ! Polar Index
real(8) :: xs ! Resultant xs
integer :: imu
real(8) :: dmu, r, f
xs = ZERO
if ((xstype == 'nu_fission' .or. xstype == 'fission' .or. xstype =='chi' &
@ -723,8 +684,6 @@ module nuclide_header
real(8) :: xs ! Resultant xs
integer :: i_azi_, i_pol_
integer :: imu
real(8) :: dmu, r, f
xs = ZERO

38
src/output.F90
View file

@ -102,10 +102,9 @@ contains
!===============================================================================
subroutine header(msg, unit, level)
character(*), intent(in) :: msg ! header message
integer, optional :: unit ! unit to write to
integer, optional :: level ! specified header level
character(*), intent(in) :: msg ! header message
integer, intent(in), optional :: unit ! unit to write to
integer, intent(in), optional :: level ! specified header level
integer :: n ! number of = signs on left
integer :: m ! number of = signs on right
@ -197,9 +196,8 @@ contains
!===============================================================================
subroutine write_message(message, level)
character(*) :: message
integer, optional :: level ! verbosity level
character(*), intent(in) :: message ! message to write
integer, intent(in), optional :: level ! verbosity level
integer :: i_start ! starting position
integer :: i_end ! ending position
@ -252,7 +250,6 @@ contains
!===============================================================================
subroutine print_particle(p)
type(Particle), intent(in) :: p
integer :: i ! index for coordinate levels
@ -486,7 +483,7 @@ contains
subroutine print_plot()
integer :: i ! loop index for plots
type(ObjectPlot), pointer :: pl => null()
type(ObjectPlot), pointer :: pl
! Display header for plotting
call header("PLOTTING SUMMARY")
@ -988,7 +985,7 @@ contains
!===============================================================================
subroutine write_surface_current(t, unit_tally)
type(TallyObject), pointer :: t
type(TallyObject), intent(in) :: t
integer, intent(in) :: unit_tally
integer :: i ! mesh index for x
@ -1160,10 +1157,9 @@ contains
!===============================================================================
function get_label(t, i_filter) result(label)
type(TallyObject), pointer :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(100) :: label ! user-specified identifier
type(TallyObject), intent(in) :: t ! tally object
integer, intent(in) :: i_filter ! index in filters array
character(100) :: label ! user-specified identifier
integer :: i ! index in cells/surfaces/etc array
integer :: bin
@ -1227,12 +1223,12 @@ contains
recursive subroutine find_offset(map, goal, univ, final, offset, path)
integer, intent(in) :: map ! Index in maps vector
integer, intent(in) :: goal ! The target cell ID
type(Universe), pointer, intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(100) :: path ! Path to offset
integer, intent(in) :: map ! Index in maps vector
integer, intent(in) :: goal ! The target cell ID
type(Universe), intent(in) :: univ ! Universe to begin search
integer, intent(in) :: final ! Target offset
integer, intent(inout) :: offset ! Current offset
character(*), intent(inout) :: path ! Path to offset
integer :: i, j ! Index over cells
integer :: k, l, m ! Indices in lattice
@ -1244,7 +1240,7 @@ contains
integer :: temp_offset ! Looped sum of offsets
logical :: this_cell = .false. ! Advance in this cell?
logical :: later_cell = .false. ! Fill cells after this one?
type(Cell), pointer:: c ! Pointer to current cell
type(Cell), pointer :: c ! Pointer to current cell
type(Universe), pointer :: next_univ ! Next universe to loop through
class(Lattice), pointer :: lat ! Pointer to current lattice

104
src/physics.F90
View file

@ -187,7 +187,6 @@ contains
!===============================================================================
subroutine sample_fission(i_nuclide, i_reaction)
integer, intent(in) :: i_nuclide ! index in nuclides array
integer, intent(out) :: i_reaction ! index in nuc % reactions array
@ -197,7 +196,6 @@ contains
real(8) :: prob
real(8) :: cutoff
type(Nuclide_CE), pointer :: nuc
type(Reaction), pointer :: rxn
! Get pointer to nuclide
nuc => nuclides(i_nuclide)
@ -222,14 +220,15 @@ contains
FISSION_REACTION_LOOP: do i = 1, nuc % n_fission
i_reaction = nuc % index_fission(i)
rxn => nuc % reactions(i_reaction)
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
associate (rxn => nuc % reactions(i_reaction))
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
end associate
! Create fission bank sites if fission occurs
if (prob > cutoff) exit FISSION_REACTION_LOOP
@ -242,9 +241,8 @@ contains
!===============================================================================
subroutine absorption(p, i_nuclide)
type(Particle), intent(inout) :: p
integer, intent(in) :: i_nuclide
integer, intent(in) :: i_nuclide
if (survival_biasing) then
! Determine weight absorbed in survival biasing
@ -283,7 +281,6 @@ contains
!===============================================================================
subroutine scatter(p, i_nuclide, i_nuc_mat)
type(Particle), intent(inout) :: p
integer, intent(in) :: i_nuclide
integer, intent(in) :: i_nuc_mat
@ -293,11 +290,10 @@ contains
real(8) :: f
real(8) :: prob
real(8) :: cutoff
type(Nuclide_CE), pointer :: nuc
type(Reaction), pointer :: rxn
real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering
real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering
real(8) :: phi ! azimuthal angle for iso-in-lab scattering
type(Nuclide_CE), pointer :: nuc
! copy incoming direction
uvw_old(:) = p % coord(1) % uvw
@ -324,11 +320,8 @@ contains
p % E, p % coord(1) % uvw, p % mu)
else
! get pointer to elastic scattering reaction
rxn => nuc % reactions(1)
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, rxn, &
call elastic_scatter(i_nuclide, nuc % reactions(1), &
p % E, p % coord(1) % uvw, p % mu, p % wgt)
end if
@ -351,28 +344,28 @@ contains
&// trim(nuc % name))
end if
rxn => nuc % reactions(i)
associate (rxn => nuc % reactions(i))
! Skip fission reactions
if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF &
.or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle
! Skip fission reactions
if (rxn % MT == N_FISSION .or. rxn % MT == N_F .or. rxn % MT == N_NF &
.or. rxn % MT == N_2NF .or. rxn % MT == N_3NF) cycle
! some materials have gas production cross sections with MT > 200 that
! are duplicates. Also MT=4 is total level inelastic scattering which
! should be skipped
if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle
! some materials have gas production cross sections with MT > 200 that
! are duplicates. Also MT=4 is total level inelastic scattering which
! should be skipped
if (rxn % MT >= 200 .or. rxn % MT == N_LEVEL) cycle
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
! if energy is below threshold for this reaction, skip it
if (i_grid < rxn % threshold) cycle
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
! add to cumulative probability
prob = prob + ((ONE - f)*rxn%sigma(i_grid - rxn%threshold + 1) &
+ f*(rxn%sigma(i_grid - rxn%threshold + 2)))
end associate
end do
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, rxn, p)
p % event_MT = rxn % MT
call inelastic_scatter(nuc, nuc%reactions(i), p)
p % event_MT = nuc%reactions(i)%MT
end if
@ -401,9 +394,8 @@ contains
!===============================================================================
subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt)
integer, intent(in) :: i_nuclide
type(Reaction), pointer :: rxn
type(Reaction), intent(in) :: rxn
real(8), intent(inout) :: E
real(8), intent(inout) :: uvw(3)
real(8), intent(out) :: mu_lab
@ -482,7 +474,6 @@ contains
!===============================================================================
subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu)
integer, intent(in) :: i_nuclide ! index in micro_xs
integer, intent(in) :: i_sab ! index in sab_tables
real(8), intent(inout) :: E ! incoming/outgoing energy
@ -740,9 +731,7 @@ contains
!===============================================================================
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff)
type(Nuclide_CE), pointer :: nuc ! target nuclide at temperature T
type(Nuclide_CE), intent(in) :: nuc ! target nuclide at temperature T
real(8), intent(out) :: v_target(3) ! target velocity
real(8), intent(in) :: v_neut(3) ! neutron velocity
real(8), intent(in) :: E ! particle energy
@ -987,11 +976,10 @@ contains
!===============================================================================
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw)
type(Nuclide_CE), pointer :: nuc ! target nuclide at temperature
real(8), intent(out) :: v_target(3)
real(8), intent(in) :: E
real(8), intent(in) :: uvw(3)
type(Nuclide_CE), intent(in) :: nuc ! target nuclide at temperature
real(8), intent(out) :: v_target(3)
real(8), intent(in) :: E
real(8), intent(in) :: uvw(3)
real(8) :: kT ! equilibrium temperature of target in MeV
real(8) :: awr ! target/neutron mass ratio
@ -1061,10 +1049,9 @@ contains
!===============================================================================
subroutine create_fission_sites(p, i_nuclide, i_reaction)
type(Particle), intent(inout) :: p
integer, intent(in) :: i_nuclide
integer, intent(in) :: i_reaction
integer, intent(in) :: i_nuclide
integer, intent(in) :: i_reaction
integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
integer :: i ! loop index
@ -1076,11 +1063,9 @@ contains
real(8) :: weight ! weight adjustment for ufs method
logical :: in_mesh ! source site in ufs mesh?
type(Nuclide_CE), pointer :: nuc
type(Reaction), pointer :: rxn
! Get pointers
nuc => nuclides(i_nuclide)
rxn => nuc % reactions(i_reaction)
! TODO: Heat generation from fission
@ -1151,7 +1136,8 @@ contains
! Sample secondary energy distribution for fission reaction and set energy
! in fission bank
fission_bank(i) % E = sample_fission_energy(nuc, rxn, p)
fission_bank(i) % E = sample_fission_energy(nuc, nuc%reactions(&
i_reaction), p)
! Set the delayed group of the neutron
fission_bank(i) % delayed_group = p % delayed_group
@ -1178,10 +1164,10 @@ contains
function sample_fission_energy(nuc, rxn, p) result(E_out)
type(Nuclide_CE), pointer :: nuc
type(Reaction), pointer :: rxn
type(Nuclide_CE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p ! Particle causing fission
real(8) :: E_out ! outgoing E of fission neutron
real(8) :: E_out ! outgoing energy of fission neutron
integer :: j ! index on nu energy grid / precursor group
integer :: lc ! index before start of energies/nu values
@ -1304,9 +1290,9 @@ contains
!===============================================================================
subroutine inelastic_scatter(nuc, rxn, p)
type(Nuclide_CE), pointer :: nuc
type(Reaction), pointer :: rxn
type(Particle), intent(inout) :: p
type(Nuclide_CE), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p
integer :: i ! loop index
integer :: law ! secondary energy distribution law
@ -1390,8 +1376,7 @@ contains
!===============================================================================
function sample_angle(rxn, E) result(mu)
type(Reaction), pointer :: rxn ! reaction
type(Reaction), intent(in) :: rxn ! reaction
real(8), intent(in) :: E ! incoming energy
real(8) :: xi ! random number on [0,1)
@ -1518,8 +1503,7 @@ contains
!===============================================================================
recursive subroutine sample_energy(edist, E_in, E_out, mu_out, A, Q)
type(DistEnergy), pointer :: edist
type(DistEnergy), intent(in) :: edist
real(8), intent(in) :: E_in ! incoming energy of neutron
real(8), intent(out) :: E_out ! outgoing energy
real(8), intent(inout), optional :: mu_out ! outgoing cosine of angle

27
src/search.F90
View file

@ -18,7 +18,7 @@ contains
! value lies in the array. This is used extensively for energy grid searching
!===============================================================================
function binary_search_real(array, n, val) result(array_index)
pure function binary_search_real(array, n, val) result(array_index)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
@ -33,7 +33,8 @@ contains
R = n
if (val < array(L) .or. val > array(R)) then
call fatal_error("Value outside of array during binary search")
array_index = -1
return
end if
n_iteration = 0
@ -49,8 +50,8 @@ contains
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
array_index = -2
return
end if
end do
@ -58,7 +59,7 @@ contains
end function binary_search_real
function binary_search_int4(array, n, val) result(array_index)
pure function binary_search_int4(array, n, val) result(array_index)
integer, intent(in) :: n
integer, intent(in) :: array(n)
@ -73,7 +74,8 @@ contains
R = n
if (val < array(L) .or. val > array(R)) then
call fatal_error("Value outside of array during binary search")
array_index = -1
return
end if
n_iteration = 0
@ -89,8 +91,8 @@ contains
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
array_index = -2
return
end if
end do
@ -98,7 +100,7 @@ contains
end function binary_search_int4
function binary_search_int8(array, n, val) result(array_index)
pure function binary_search_int8(array, n, val) result(array_index)
integer, intent(in) :: n
integer(8), intent(in) :: array(n)
@ -113,7 +115,8 @@ contains
R = n
if (val < array(L) .or. val > array(R)) then
call fatal_error("Value outside of array during binary search")
array_index = -1
return
end if
n_iteration = 0
@ -129,8 +132,8 @@ contains
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
call fatal_error("Reached maximum number of iterations on binary &
&search.")
array_index = -2
return
end if
end do

31
src/simple_string.F90
View file

@ -14,8 +14,7 @@ contains
! TO_LOWER converts a string to all lower case characters
!===============================================================================
function to_lower(word) result(word_lower)
pure function to_lower(word) result(word_lower)
character(*), intent(in) :: word
character(len=len(word)) :: word_lower
@ -37,8 +36,7 @@ contains
! TO_UPPER converts a string to all upper case characters
!===============================================================================
function to_upper(word) result(word_upper)
pure function to_upper(word) result(word_upper)
character(*), intent(in) :: word
character(len=len(word)) :: word_upper
@ -60,8 +58,7 @@ contains
! IS_NUMBER determines whether a string of characters is all 0-9 characters
!===============================================================================
function is_number(word) result(number)
pure function is_number(word) result(number)
character(*), intent(in) :: word
logical :: number
@ -81,10 +78,9 @@ contains
! sequence of characters
!===============================================================================
logical function starts_with(str, seq)
character(*) :: str ! string to check
character(*) :: seq ! sequence of characters
pure logical function starts_with(str, seq)
character(*), intent(in) :: str ! string to check
character(*), intent(in) :: seq ! sequence of characters
integer :: i
integer :: i_start
@ -116,10 +112,9 @@ contains
! of characters
!===============================================================================
logical function ends_with(str, seq)
character(*) :: str ! string to check
character(*) :: seq ! sequence of characters
pure logical function ends_with(str, seq)
character(*), intent(in) :: str ! string to check
character(*), intent(in) :: seq ! sequence of characters
integer :: i_start
integer :: str_len
@ -145,7 +140,7 @@ contains
! integer.
!===============================================================================
function count_digits(num) result(n_digits)
pure function count_digits(num) result(n_digits)
integer, intent(in) :: num
integer :: n_digits
@ -163,7 +158,7 @@ contains
! INT4_TO_STR converts an integer(4) to a string.
!===============================================================================
function int4_to_str(num) result(str)
pure function int4_to_str(num) result(str)
integer, intent(in) :: num
character(11) :: str
@ -177,7 +172,7 @@ contains
! INT8_TO_STR converts an integer(8) to a string.
!===============================================================================
function int8_to_str(num) result(str)
pure function int8_to_str(num) result(str)
integer(8), intent(in) :: num
character(21) :: str
@ -193,7 +188,7 @@ contains
! are used.
!===============================================================================
function real_to_str(num, sig_digits) result(string)
pure function real_to_str(num, sig_digits) result(string)
real(8), intent(in) :: num ! number to convert
integer, optional, intent(in) :: sig_digits ! # of significant digits

3
src/source.F90
View file

@ -97,8 +97,7 @@ contains
!===============================================================================
subroutine sample_external_source(site)
type(Bank), pointer :: site ! source site
type(Bank), intent(inout) :: site ! source site
integer :: i ! dummy loop index
real(8) :: r(3) ! sampled coordinates

7
src/string.F90
View file

@ -22,7 +22,6 @@ contains
!===============================================================================
subroutine split_string(string, words, n)
character(*), intent(in) :: string
character(*), intent(out) :: words(MAX_WORDS)
integer, intent(out) :: n
@ -163,7 +162,7 @@ contains
! string = concatenated string
!===============================================================================
function concatenate(words, n_words) result(string)
pure function concatenate(words, n_words) result(string)
integer, intent(in) :: n_words
character(*), intent(in) :: words(n_words)
@ -218,7 +217,7 @@ contains
! STR_TO_INT converts a string to an integer.
!===============================================================================
function str_to_int(str) result(num)
pure function str_to_int(str) result(num)
character(*), intent(in) :: str
integer(8) :: num
@ -243,7 +242,7 @@ contains
! STR_TO_REAL converts an arbitrary string to a real(8)
!===============================================================================
function str_to_real(string) result(num)
pure function str_to_real(string) result(num)
character(*), intent(in) :: string
real(8) :: num

489
src/tally.F90
View file

@ -5,7 +5,7 @@ module tally
use error, only: fatal_error
use geometry_header
use global
use math, only: t_percentile, calc_pn, calc_rn, evaluate_legendre
use math, only: t_percentile, calc_pn, calc_rn
use mesh, only: get_mesh_bin, bin_to_mesh_indices, &
get_mesh_indices, mesh_indices_to_bin, &
mesh_intersects_2d, mesh_intersects_3d
@ -34,20 +34,20 @@ module tally
atom_density, flux)
import Particle
import TallyObject
type(Particle), intent(in) :: p
type(TallyObject), pointer, intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
end subroutine score_general_intfc
subroutine get_scoring_bins_intfc(p, i_tally, found_bin)
import Particle
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
end subroutine get_scoring_bins_intfc
end interface
@ -79,13 +79,13 @@ contains
subroutine score_general_ce(p, t, start_index, filter_index, i_nuclide, &
atom_density, flux)
type(Particle), intent(in) :: p
type(TallyObject), pointer, intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
integer :: i ! loop index for scoring bins
integer :: l ! loop index for nuclides in material
@ -104,11 +104,7 @@ contains
real(8) :: score ! analog tally score
real(8) :: macro_total ! material macro total xs
real(8) :: macro_scatt ! material macro scatt xs
real(8) :: uvw(3) ! particle direction
real(8) :: E ! particle energy
type(Material), pointer :: mat
type(Reaction), pointer :: rxn
type(Nuclide_CE), pointer :: nuc
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
@ -259,24 +255,20 @@ contains
! of one.
score = p % last_wgt
else
do m = 1, nuclides(p % event_nuclide) % n_reaction
! Check if this is the desired MT
if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then
! Found the reaction, set our pointer and move on with life
rxn => nuclides(p % event_nuclide) % reactions(m)
exit
end if
end do
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end associate
end if
@ -296,24 +288,20 @@ contains
! of one.
score = p % last_wgt
else
do m = 1, nuclides(p % event_nuclide) % n_reaction
! Check if this is the desired MT
if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then
! Found the reaction, set our pointer and move on with life
rxn => nuclides(p % event_nuclide) % reactions(m)
exit
end if
end do
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end associate
end if
@ -333,24 +321,20 @@ contains
! of one.
score = p % last_wgt
else
do m = 1, nuclides(p % event_nuclide) % n_reaction
! Check if this is the desired MT
if (p % event_MT == nuclides(p % event_nuclide) % reactions(m) % MT) then
! Found the reaction, set our pointer and move on with life
rxn => nuclides(p % event_nuclide) % reactions(m)
exit
end if
end do
m = nuclides(p%event_nuclide)%reaction_index% &
get_key(p % event_MT)
! Get multiplicity and apply to score
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
associate (rxn => nuclides(p%event_nuclide)%reactions(m))
if (rxn % multiplicity_with_E) then
! Then the multiplicity was already incorporated in to p % wgt
! per the scattering routine,
score = p % wgt
else
! Grab the multiplicity from the rxn
score = p % last_wgt * rxn % multiplicity
end if
end associate
end if
@ -505,9 +489,6 @@ contains
! delayed-nu-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
! Get the event nuclide
nuc => nuclides(p % event_nuclide)
! Check if the delayed group filter is present
if (dg_filter > 0) then
@ -519,11 +500,11 @@ contains
d = t % filters(dg_filter) % int_bins(d_bin)
! Compute the yield for this delayed group
yield = yield_delayed(nuc, E, d)
yield = yield_delayed(nuclides(p % event_nuclide), E, d)
! Compute the score and tally to bin
score = p % absorb_wgt * yield * micro_xs(p % event_nuclide) &
% fission * nu_delayed(nuc, E) / &
% fission * nu_delayed(nuclides(p % event_nuclide), E) / &
micro_xs(p % event_nuclide) % absorption
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
@ -533,7 +514,7 @@ contains
! by multiplying the absorbed weight by the fraction of the
! delayed-nu-fission xs to the absorption xs
score = p % absorb_wgt * micro_xs(p % event_nuclide) &
% fission * nu_delayed(nuc, E) / &
% fission * nu_delayed(nuclides(p % event_nuclide), E) / &
micro_xs(p % event_nuclide) % absorption
end if
end if
@ -574,9 +555,6 @@ contains
! Check if tally is on a single nuclide
if (i_nuclide > 0) then
! Get the nuclide of interest
nuc => nuclides(i_nuclide)
! Check if the delayed group filter is present
if (dg_filter > 0) then
@ -587,11 +565,11 @@ contains
d = t % filters(dg_filter) % int_bins(d_bin)
! Compute the yield for this delayed group
yield = yield_delayed(nuc, E, d)
yield = yield_delayed(nuclides(i_nuclide), E, d)
! Compute the score and tally to bin
score = micro_xs(i_nuclide) % fission * yield &
* nu_delayed(nuc, E) * atom_density * flux
* nu_delayed(nuclides(i_nuclide), E) * atom_density * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
cycle SCORE_LOOP
@ -599,27 +577,24 @@ contains
! If the delayed group filter is not present, compute the score
! by multiplying the delayed-nu-fission macro xs by the flux
score = micro_xs(i_nuclide) % fission * nu_delayed(nuc, E)&
* atom_density * flux
score = micro_xs(i_nuclide) % fission * &
nu_delayed(nuclides(i_nuclide), E) * atom_density * flux
end if
! Tally is on total nuclides
else
! Get pointer to current material
mat => materials(p % material)
! Check if the delayed group filter is present
if (dg_filter > 0) then
! Loop over all nuclides in the current material
do l = 1, mat % n_nuclides
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = mat % atom_density(l)
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = mat % nuclide(l)
i_nuc = materials(p % material) % nuclide(l)
! Loop over all delayed group bins and tally to them individually
do d_bin = 1, t % filters(dg_filter) % n_bins
@ -627,15 +602,12 @@ contains
! Get the delayed group for this bin
d = t % filters(dg_filter) % int_bins(d_bin)
! Get the current nuclide
nuc => nuclides(i_nuc)
! Get the yield for the desired nuclide and delayed group
yield = yield_delayed(nuc, E, d)
yield = yield_delayed(nuclides(i_nuc), E, d)
! Compute the score and tally to bin
score = micro_xs(i_nuc) % fission * yield &
* nu_delayed(nuc, E) * atom_density_ * flux
* nu_delayed(nuclides(i_nuc), E) * atom_density_ * flux
call score_fission_delayed_dg(t, d_bin, score, score_index)
end do
end do
@ -645,13 +617,13 @@ contains
score = ZERO
! Loop over all nuclides in the current material
do l = 1, mat % n_nuclides
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = mat % atom_density(l)
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = mat % nuclide(l)
i_nuc = materials(p % material) % nuclide(l)
! Accumulate the contribution from each nuclide
score = score + micro_xs(i_nuc) % fission &
@ -663,38 +635,67 @@ contains
case (SCORE_KAPPA_FISSION)
! Determine kappa-fission cross section on the fly. The ENDF standard
! (ENDF-102) states that MT 18 stores the fission energy as the Q_value
! (fission(1))
score = ZERO
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scale by kappa-fission
if (micro_xs(p % event_nuclide) % absorption > ZERO) then
score = p % absorb_wgt * &
micro_xs(p % event_nuclide) % kappa_fission / &
micro_xs(p % event_nuclide) % absorption
else
score = ZERO
end if
associate (nuc => nuclides(p%event_nuclide))
if (micro_xs(p%event_nuclide)%absorption > ZERO .and. &
nuc%fissionable) then
score = p%absorb_wgt * &
nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(p%event_nuclide)%fission / &
micro_xs(p%event_nuclide)%absorption
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
score = p % last_wgt * &
micro_xs(p % event_nuclide) % kappa_fission / &
micro_xs(p % event_nuclide) % absorption
associate (nuc => nuclides(p%event_nuclide))
if (nuc%fissionable) then
score = p%last_wgt * &
nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(p%event_nuclide)%fission / &
micro_xs(p%event_nuclide)%absorption
end if
end associate
end if
else
if (i_nuclide > 0) then
score = micro_xs(i_nuclide) % kappa_fission * atom_density * flux
associate (nuc => nuclides(i_nuclide))
if (nuc%fissionable) then
score = nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(i_nuclide)%fission * atom_density * flux
end if
end associate
else
score = material_xs % kappa_fission * flux
do l = 1, materials(p%material)%n_nuclides
! Determine atom density and index of nuclide
atom_density_ = materials(p%material)%atom_density(l)
i_nuc = materials(p%material)%nuclide(l)
! If nuclide is fissionable, accumulate kappa fission
associate(nuc => nuclides(i_nuc))
if (nuc % fissionable) then
score = score + nuc%reactions(nuc%index_fission(1))%Q_value * &
micro_xs(i_nuc)%fission * atom_density_ * flux
end if
end associate
end do
end if
end if
case (SCORE_EVENTS)
! Simply count number of scoring events
score = ONE
@ -730,14 +731,10 @@ contains
score = ZERO
if (i_nuclide > 0) then
! TODO: The following search for the matching reaction could
! be replaced by adding a dictionary on each Nuclide instance
! of the form {MT: i_reaction, ...}
REACTION_LOOP: do m = 1, nuclides(i_nuclide) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuclide) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
if (nuclides(i_nuclide)%reaction_index%has_key(score_bin)) then
m = nuclides(i_nuclide)%reaction_index%get_key(score_bin)
associate (rxn => nuclides(i_nuclide) % reactions(m))
! Retrieve index on nuclide energy grid and interpolation
! factor
i_energy = micro_xs(i_nuclide) % index_grid
@ -747,26 +744,20 @@ contains
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density * flux
end if
exit REACTION_LOOP
end if
end do REACTION_LOOP
end associate
end if
else
! Get pointer to current material
mat => materials(p % material)
do l = 1, mat % n_nuclides
do l = 1, materials(p % material) % n_nuclides
! Get atom density
atom_density_ = mat % atom_density(l)
atom_density_ = materials(p % material) % atom_density(l)
! Get index in nuclides array
i_nuc = mat % nuclide(l)
! TODO: The following search for the matching reaction could
! be replaced by adding a dictionary on each Nuclide
! instance of the form {MT: i_reaction, ...}
do m = 1, nuclides(i_nuc) % n_reaction
! Get pointer to reaction
rxn => nuclides(i_nuc) % reactions(m)
! Check if this is the desired MT
if (score_bin == rxn % MT) then
i_nuc = materials(p % material) % nuclide(l)
if (nuclides(i_nuc)%reaction_index%has_key(score_bin)) then
m = nuclides(i_nuc)%reaction_index%get_key(score_bin)
associate (rxn => nuclides(i_nuc) % reactions(m))
! Retrieve index on nuclide energy grid and interpolation
! factor
i_energy = micro_xs(i_nuc) % index_grid
@ -776,9 +767,8 @@ contains
rxn%threshold + 1) + f * rxn % sigma(i_energy - &
rxn%threshold + 2)) * atom_density_ * flux
end if
exit
end if
end do
end associate
end if
end do
end if
@ -801,26 +791,22 @@ contains
subroutine score_general_mg(p, t, start_index, filter_index, i_nuclide, &
atom_density, flux)
type(Particle), intent(in) :: p
type(TallyObject), pointer, intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: start_index
integer, intent(in) :: i_nuclide
integer, intent(in) :: filter_index ! for % results
real(8), intent(in) :: flux ! flux estimate
real(8), intent(in) :: atom_density ! atom/b-cm
integer :: i ! loop index for scoring bins
integer :: q ! loop index for scoring bins
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
integer :: score_index ! scoring bin index
integer :: g ! Group of interest
real(8) :: score ! analog tally score
real(8) :: macro_total ! material macro total xs
real(8) :: macro_scatt ! material macro scatt xs
real(8) :: micro_abs ! nuclidic microscopic abs
class(Nuclide_MG), pointer :: nuc
if (i_nuclide > 0) nuc => nuclides_MG(i_nuclide) % obj
i = 0
SCORE_LOOP: do q = 1, t % n_user_score_bins
@ -873,8 +859,10 @@ contains
else
if (i_nuclide > 0) then
score = nuc % get_xs(p % g, 'total', UVW=p % coord(i) % uvw) * &
atom_density * flux
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = nuc % get_xs(p % g, 'total', UVW=p % coord(i) % uvw) * &
atom_density * flux
end associate
else
score = material_xs % total * flux
end if
@ -913,8 +901,10 @@ contains
else
! Note SCORE_SCATTER_N not available for tracklength/collision.
if (i_nuclide > 0) then
score = nuc % get_xs(p % g, 'scatter', UVW=p % coord(i) % uvw) * &
atom_density * flux
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = nuc % get_xs(p % g, 'scatter', UVW=p % coord(i) % uvw) * &
atom_density * flux
end associate
else
! Get the scattering x/s (stored in % elastic)
score = material_xs % elastic * flux
@ -922,8 +912,10 @@ contains
end if
if (i_nuclide > 0) then
score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, &
p % last_uvw, p % mu)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, &
p % last_uvw, p % mu)
end associate
else
score = score / &
macro_xs(p % material) % obj % get_xs(p % g, 'mult', &
@ -945,8 +937,10 @@ contains
score = p % last_wgt
if (i_nuclide > 0) then
score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, &
p % last_uvw, p % mu)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, &
p % last_uvw, p % mu)
end associate
else
score = score / &
macro_xs(p % material) % obj % get_xs(p % g, 'mult', &
@ -968,8 +962,10 @@ contains
score = p % last_wgt
if (i_nuclide > 0) then
score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, &
p % last_uvw, p % mu)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = score * nuc % get_xs(p % g, 'f_mu/mult', p % last_g, &
p % last_uvw, p % mu)
end associate
else
score = score / &
macro_xs(p % material) % obj % get_xs(p % g, 'mult', &
@ -987,8 +983,10 @@ contains
! neutrons exiting a reaction with neutrons in the exit channel
score = p % wgt
if (i_nuclide > 0) then
score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, &
p % last_uvw, p % mu)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, &
p % last_uvw, p % mu)
end associate
end if
@ -1004,8 +1002,10 @@ contains
! neutrons exiting a reaction with neutrons in the exit channel
score = p % wgt
if (i_nuclide > 0) then
score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, &
p % last_uvw, p % mu)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, &
p % last_uvw, p % mu)
end associate
end if
@ -1021,8 +1021,10 @@ contains
! neutrons exiting a reaction with neutrons in the exit channel
score = p % wgt
if (i_nuclide > 0) then
score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, &
p % last_uvw, p % mu)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = score * nuc % get_xs(p % g, 'f_mu', p % last_g, &
p % last_uvw, p % mu)
end associate
end if
@ -1066,8 +1068,10 @@ contains
else
if (i_nuclide > 0) then
score = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw) &
* atom_density * flux
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw) &
* atom_density * flux
end associate
else
score = material_xs % absorption * flux
end if
@ -1080,29 +1084,35 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission
micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
if (micro_abs > ZERO) then
score = p % absorb_wgt * &
nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) &
/ micro_abs
else
score = ZERO
end if
associate (nuc => nuclides_MG(i_nuclide) % obj)
micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
if (micro_abs > ZERO) then
score = p % absorb_wgt * &
nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) &
/ micro_abs
else
score = ZERO
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for the
! fission reaction rate
score = p % last_wgt &
* nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) &
/ nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = p % last_wgt &
* nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) &
/ nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
end associate
end if
else
if (i_nuclide > 0) then
score = nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) * &
atom_density * flux
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) * &
atom_density * flux
end associate
else
score = material_xs % fission * flux
end if
@ -1126,14 +1136,16 @@ contains
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! nu-fission
micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
if (micro_abs > ZERO) then
score = p % absorb_wgt * &
nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) / &
micro_abs
else
score = ZERO
end if
associate (nuc => nuclides_MG(i_nuclide) % obj)
micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
if (micro_abs > ZERO) then
score = p % absorb_wgt * &
nuc % get_xs(p % g, 'fission', UVW=p % coord(i) % uvw) / &
micro_abs
else
score = ZERO
end if
end associate
else
! Skip any non-fission events
if (.not. p % fission) cycle SCORE_LOOP
@ -1147,8 +1159,10 @@ contains
else
if (i_nuclide > 0) then
score = nuc % get_xs(p % g, 'nu_fission', UVW=p % coord(i) % uvw) &
* atom_density * flux
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = nuc % get_xs(p % g, 'nu_fission', UVW=p % coord(i) % uvw) &
* atom_density * flux
end associate
else
score = material_xs % nu_fission * flux
end if
@ -1156,36 +1170,43 @@ contains
case (SCORE_KAPPA_FISSION)
! Determine kappa-fission cross section
score = ZERO
if (t % estimator == ESTIMATOR_ANALOG) then
if (survival_biasing) then
! No fission events occur if survival biasing is on -- need to
! calculate fraction of absorptions that would have resulted in
! fission scale by kappa-fission
micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
if (micro_abs > ZERO) then
score = p % absorb_wgt * &
nuc % get_xs(p % g, 'k_fission', UVW=p % coord(i) % uvw) / &
micro_abs
else
score = ZERO
end if
associate (nuc => nuclides_MG(i_nuclide) % obj)
micro_abs = nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
if (micro_abs > ZERO) then
score = p % absorb_wgt * &
nuc % get_xs(p % g, 'k_fission', UVW=p % coord(i) % uvw) / &
micro_abs
end if
end associate
else
! Skip any non-absorption events
if (p % event == EVENT_SCATTER) cycle SCORE_LOOP
! All fission events will contribute, so again we can use
! particle's weight entering the collision as the estimate for
! the fission energy production rate
score = p % last_wgt * &
nuc % get_xs(p % g, 'k_fission', UVW=p % coord(i) % uvw) / &
nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = p % last_wgt * &
nuc % get_xs(p % g, 'k_fission', UVW=p % coord(i) % uvw) / &
nuc % get_xs(p % g, 'absorption', UVW=p % coord(i) % uvw)
end associate
end if
else
if (i_nuclide > 0) then
score = nuc % get_xs(p % g, 'k_fission', UVW=p % coord(i) % uvw) &
* atom_density * flux
associate (nuc => nuclides_MG(i_nuclide) % obj)
score = nuc % get_xs(p % g, 'k_fission', UVW=p % coord(i) % uvw) &
* atom_density * flux
end associate
else
score = material_xs % kappa_fission * flux
score = macro_xs(p % material) % obj % get_xs(p % g, 'k_fission', &
UVW=p % coord(i) % uvw)
end if
end if
@ -1212,13 +1233,13 @@ contains
subroutine expand_and_score(p, t, score_index, filter_index, score_bin, &
score, i)
type(Particle), intent(in) :: p
type(TallyObject), pointer, intent(inout) :: t
integer, intent(inout) :: score_index
integer, intent(in) :: filter_index ! for % results
integer, intent(in) :: score_bin ! score of concern
real(8), intent(inout) :: score ! data to score
integer, intent(inout) :: i ! Working index
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(inout) :: score_index
integer, intent(in) :: filter_index ! for % results
integer, intent(in) :: score_bin ! score of concern
real(8), intent(inout) :: score ! data to score
integer, intent(inout) :: i ! Working index
integer :: num_nm ! Number of N,M orders in harmonic
integer :: n ! Moment loop index
@ -1327,9 +1348,9 @@ contains
subroutine score_all_nuclides(p, i_tally, flux, filter_index)
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
real(8), intent(in) :: flux
integer, intent(in) :: filter_index
integer, intent(in) :: i_tally
real(8), intent(in) :: flux
integer, intent(in) :: filter_index
integer :: i ! loop index for nuclides in material
integer :: i_nuclide ! index in nuclides array
@ -1487,10 +1508,9 @@ contains
!===============================================================================
subroutine score_fission_eout_ce(p, t, i_score)
type(Particle), intent(in) :: p
type(TallyObject), pointer :: t
integer, intent(in) :: i_score ! index for score
type(TallyObject), intent(inout) :: t
integer, intent(in) :: i_score ! index for score
integer :: i ! index of outgoing energy filter
integer :: n ! number of energies on filter
@ -1541,11 +1561,10 @@ contains
end subroutine score_fission_eout_ce
subroutine score_fission_eout_mg(p, t, i_score)
type(Particle), intent(in) :: p
type(TallyObject), pointer :: t
integer, intent(in) :: i_score ! index for score
subroutine score_fission_eout_mg(p, t, i_score)
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: i_score ! index for score
integer :: i ! index of outgoing energy filter
integer :: n ! number of energies on filter
@ -1614,7 +1633,7 @@ contains
subroutine score_fission_delayed_eout(p, t, i_score)
type(Particle), intent(in) :: p
type(Particle), intent(in) :: p
type(TallyObject), intent(inout) :: t
integer, intent(in) :: i_score ! index for score
@ -1741,7 +1760,7 @@ contains
subroutine score_tracklength_tally(p, distance)
type(Particle), intent(in) :: p
real(8), intent(in) :: distance
real(8), intent(in) :: distance
integer :: i
integer :: i_tally
@ -1854,8 +1873,8 @@ contains
subroutine score_tl_on_mesh(p, i_tally, d_track)
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
real(8), intent(in) :: d_track
integer, intent(in) :: i_tally
real(8), intent(in) :: d_track
integer :: i ! loop index for filter/score bins
integer :: j ! loop index for direction
@ -1881,9 +1900,9 @@ contains
logical :: end_in_mesh ! ending coordinates inside mesh?
real(8) :: theta
real(8) :: phi
type(TallyObject), pointer :: t
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
type(Material), pointer :: mat
type(Material), pointer :: mat
t => tallies(i_tally)
matching_bins(1:t%n_filters) = 1
@ -2247,9 +2266,9 @@ contains
subroutine get_scoring_bins_ce(p, i_tally, found_bin)
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
integer :: i ! loop index for filters
integer :: j
@ -2257,7 +2276,7 @@ contains
integer :: offset ! offset for distribcell
real(8) :: E ! particle energy
real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively
type(TallyObject), pointer :: t
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
found_bin = .true.
@ -2452,9 +2471,9 @@ contains
subroutine get_scoring_bins_mg(p, i_tally, found_bin)
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
type(Particle), intent(in) :: p
integer, intent(in) :: i_tally
logical, intent(out) :: found_bin
integer :: i ! loop index for filters
integer :: j
@ -2694,7 +2713,7 @@ contains
logical :: x_same ! same starting/ending x index (i)
logical :: y_same ! same starting/ending y index (j)
logical :: z_same ! same starting/ending z index (k)
type(TallyObject), pointer :: t
type(TallyObject), pointer :: t
type(RegularMesh), pointer :: m
TALLY_LOOP: do i = 1, active_current_tallies % size()

79
src/tally_header.F90
View file

@ -58,10 +58,6 @@ module tally_header
integer :: offset = 0 ! Only used for distribcell filters
integer, allocatable :: int_bins(:)
real(8), allocatable :: real_bins(:) ! Only used for energy filters
! Type-Bound procedures
contains
procedure :: clear => tallyfilter_clear ! Deallocates TallyFilter
end type TallyFilter
!===============================================================================
@ -133,81 +129,6 @@ module tally_header
! Multi-Group Specific Information To Enable Rapid Tallying
logical :: energy_matches_groups = .false.
logical :: energyout_matches_groups = .false.
! Type-Bound procedures
contains
procedure :: clear => tallyobject_clear ! Deallocates TallyObject
end type TallyObject
contains
!===============================================================================
! TALLYFILTER_CLEAR deallocates a TallyFilter element and sets it to its as
! initialized state.
!===============================================================================
subroutine tallyfilter_clear(this)
class(TallyFilter), intent(inout) :: this ! The TallyFilter to be cleared
this % type = NONE
this % n_bins = 0
if (allocated(this % int_bins)) &
deallocate(this % int_bins)
if (allocated(this % real_bins)) &
deallocate(this % real_bins)
end subroutine tallyfilter_clear
!===============================================================================
! TALLYOBJECT_CLEAR deallocates a TallyObject element and sets it to its as
! initialized state.
!===============================================================================
subroutine tallyobject_clear(this)
class(TallyObject), intent(inout) :: this ! The TallyObject to be cleared
integer :: i ! Loop Index
! This routine will go through each item in TallyObject and set the value
! to its default, as-initialized values, including deallocations.
this % name = ""
if (allocated(this % filters)) then
do i = 1, size(this % filters)
call this % filters(i) % clear()
end do
deallocate(this % filters)
end if
if (allocated(this % stride)) &
deallocate(this % stride)
this % find_filter = 0
this % n_nuclide_bins = 0
if (allocated(this % nuclide_bins)) &
deallocate(this % nuclide_bins)
this % all_nuclides = .false.
this % n_score_bins = 0
if (allocated(this % score_bins)) &
deallocate(this % score_bins)
if (allocated(this % moment_order)) &
deallocate(this % moment_order)
this % n_user_score_bins = 0
if (allocated(this % results)) &
deallocate(this % results)
this % reset = .false.
this % n_realizations = 0
if (allocated(this % triggers)) &
deallocate (this % triggers)
this % n_triggers = 0
end subroutine tallyobject_clear
end module tally_header

4
src/tally_initialize.F90
View file

@ -38,7 +38,7 @@ contains
integer :: j ! loop index for filters
integer :: n ! temporary stride
integer :: max_n_filters = 0 ! maximum number of filters
type(TallyObject), pointer :: t => null()
type(TallyObject), pointer :: t
TALLY_LOOP: do i = 1, n_tallies
! Get pointer to tally
@ -88,7 +88,7 @@ contains
integer :: k ! loop index for bins
integer :: bin ! filter bin entries
integer :: type ! type of tally filter
type(TallyObject), pointer :: t => null()
type(TallyObject), pointer :: t
! allocate tally map array -- note that we don't need a tally map for the
! energy_in and energy_out filters

12
src/timer_header.F90
View file

@ -29,13 +29,11 @@ contains
!===============================================================================
subroutine timer_start(self)
class(Timer), intent(inout) :: self
! Turn timer on and measure starting time
self % running = .true.
call system_clock(self % start_counts)
end subroutine timer_start
!===============================================================================
@ -43,7 +41,6 @@ contains
!===============================================================================
function timer_get_value(self) result(elapsed)
class(Timer), intent(in) :: self ! the timer
real(8) :: elapsed ! total elapsed time
@ -58,7 +55,6 @@ contains
else
elapsed = self % elapsed
end if
end function timer_get_value
!===============================================================================
@ -66,30 +62,26 @@ contains
!===============================================================================
subroutine timer_stop(self)
class(Timer), intent(inout) :: self
! Check to make sure timer was running
if (.not. self % running) return
! Stop timer and add time
self % elapsed = timer_get_value(self)
self % elapsed = self % get_value()
self % running = .false.
end subroutine timer_stop
!===============================================================================
! TIMER_RESET resets a timer to have a zero value
!===============================================================================
subroutine timer_reset(self)
pure subroutine timer_reset(self)
class(Timer), intent(inout) :: self
self % running = .false.
self % start_counts = 0
self % elapsed = ZERO
end subroutine timer_reset
end module timer_header

10
src/trigger_header.F90
View file

@ -1,6 +1,6 @@
module trigger_header
use constants, only: NONE, N_FILTER_TYPES
use constants, only: NONE, N_FILTER_TYPES, ZERO
implicit none
@ -13,9 +13,9 @@ module trigger_header
real(8) :: threshold ! a convergence threshold
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) :: rel_err =0.0 ! temp rel. err. container
real(8) :: variance = ZERO ! temp variance container
real(8) :: std_dev = ZERO ! temp std. dev. container
real(8) :: rel_err = ZERO ! temp rel. err. container
end type TriggerObject
!===============================================================================
@ -23,7 +23,7 @@ module trigger_header
!===============================================================================
type KTrigger
integer :: trigger_type = 0
real(8) :: threshold = 0
real(8) :: threshold = ZERO
end type KTrigger
end module trigger_header

2
src/xml_interface.F90
View file

@ -117,7 +117,7 @@ contains
type(Node), pointer, intent(out) :: out_ptr
logical :: found_
type(NodeList), pointer :: elem_list => null()
type(NodeList), pointer :: elem_list
! Set found to false
found_ = .false.

4
tests/test_many_scores/results_true.dat
View file

@ -33,8 +33,8 @@ tally 1:
7.620560E-01
1.816851E+00
1.102658E+00
1.338067E+02
5.986137E+03
1.337996E+02
5.985519E+03
2.247257E+01
1.683779E+02
1.512960E-01

20
tests/test_score_kappafission/results_true.dat
View file

@ -1,17 +1,17 @@
k-combined:
9.903196E-01 4.279617E-02
tally 1:
2.266048E+02
1.049833E+04
2.266169E+02
1.049923E+04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.366139E+02
3.859561E+03
1.366590E+02
3.861651E+03
tally 2:
2.402814E+02
1.174003E+04
2.403775E+02
1.175130E+04
0.000000E+00
0.000000E+00
0.000000E+00
@ -19,11 +19,11 @@ tally 2:
1.270420E+02
3.297537E+03
tally 3:
2.217075E+02
1.003168E+04
2.217588E+02
1.003581E+04
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.375693E+02
3.872389E+03
1.376303E+02
3.875598E+03